Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2018-2019 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 hydro_energy
21 :
22 : use star_private_def
23 : use const_def, only: dp, ln10, pi, pi4
24 : use utils_lib, only: mesa_error, is_bad
25 : use auto_diff
26 : use auto_diff_support
27 : use star_utils, only: em1, e00, ep1, set_energy_eqn_scal
28 :
29 : implicit none
30 :
31 : private
32 : public :: do1_energy_eqn
33 :
34 : contains
35 :
36 52348 : subroutine do1_energy_eqn( & ! energy conservation
37 : s, k, do_chem, nvar, ierr)
38 : use star_utils, only: store_partials
39 : type (star_info), pointer :: s
40 : integer, intent(in) :: k, nvar
41 : logical, intent(in) :: do_chem
42 : integer, intent(out) :: ierr
43 1936876 : real(dp), dimension(nvar) :: d_dm1, d_d00, d_dp1
44 : include 'formats'
45 : call get1_energy_eqn( &
46 : s, k, do_chem, nvar, &
47 52348 : d_dm1, d_d00, d_dp1, ierr)
48 52348 : if (ierr /= 0) then
49 0 : if (s% report_ierr) write(*,2) 'ierr /= 0 for get1_energy_eqn', k
50 : return
51 : end if
52 : call store_partials( &
53 52348 : s, k, s% i_dlnE_dt, nvar, d_dm1, d_d00, d_dp1, 'do1_energy_eqn', ierr)
54 52348 : end subroutine do1_energy_eqn
55 :
56 :
57 52348 : subroutine get1_energy_eqn( &
58 52348 : s, k, do_chem, nvar, d_dm1, d_d00, d_dp1, ierr)
59 52348 : use eos_def, only: i_grad_ad, i_lnPgas, i_lnE
60 : use eps_grav, only: eval_eps_grav_and_partials
61 : use accurate_sum_auto_diff_star_order1
62 : use auto_diff_support
63 : type (star_info), pointer :: s
64 : integer, intent(in) :: k, nvar
65 : logical, intent(in) :: do_chem
66 : real(dp), intent(out), dimension(nvar) :: d_dm1, d_d00, d_dp1
67 : integer, intent(out) :: ierr
68 :
69 : type(auto_diff_real_star_order1) :: resid_ad, &
70 : dL_dm_ad, sources_ad, others_ad, d_turbulent_energy_dt_ad, &
71 : dwork_dm_ad, eps_grav_ad, dke_dt_ad, dpe_dt_ad, de_dt_ad
72 : type(accurate_auto_diff_real_star_order1) :: esum_ad
73 52348 : real(dp) :: residual, dm, dt, scal
74 : real(dp), dimension(s% species) :: &
75 1308700 : d_dwork_dxam1, d_dwork_dxa00, d_dwork_dxap1
76 : integer :: nz, i_dlnE_dt, i_lum, i_v
77 : logical :: test_partials, doing_op_split_burn, eps_grav_form
78 :
79 : include 'formats'
80 :
81 : !test_partials = (k == s% solver_test_partials_k)
82 52348 : test_partials = .false.
83 :
84 52348 : ierr = 0
85 52348 : call init
86 :
87 52348 : call setup_eps_grav(ierr); if (ierr /= 0) return ! do this first - it sets eps_grav_form
88 52348 : call setup_de_dt_and_friends(ierr); if (ierr /= 0) return
89 52348 : call setup_dwork_dm(ierr); if (ierr /= 0) return
90 52348 : call setup_dL_dm(ierr); if (ierr /= 0) return
91 52348 : call setup_sources_and_others(ierr); if (ierr /= 0) return
92 52348 : call setup_d_turbulent_energy_dt(ierr); if (ierr /= 0) return
93 52348 : call set_energy_eqn_scal(s, k, scal, ierr); if (ierr /= 0) return
94 :
95 52348 : s% dL_dm(k) = dL_dm_ad%val
96 52348 : s% dwork_dm(k) = dwork_dm_ad%val
97 52348 : s% energy_sources(k) = sources_ad%val
98 : ! nuclear heating, non_nuc_neu_cooling, irradiation heating, extra_heat, eps_mdot
99 52348 : s% energy_others(k) = others_ad%val
100 : ! eps_WD_sedimentation, eps_diffusion, eps_pre_mix, eps_phase_separation
101 : ! sum terms in esum_ad using accurate_auto_diff_real_star_order1
102 52348 : if (eps_grav_form) then ! for this case, dwork_dm doesn't include work by P since that is in eps_grav
103 : esum_ad = - dL_dm_ad + sources_ad + &
104 0 : others_ad - d_turbulent_energy_dt_ad - dwork_dm_ad + eps_grav_ad
105 52348 : else if (s% using_velocity_time_centering .and. &
106 : s% use_P_d_1_div_rho_form_of_work_when_time_centering_velocity) then
107 : esum_ad = - dL_dm_ad + sources_ad + &
108 0 : others_ad - d_turbulent_energy_dt_ad - dwork_dm_ad - de_dt_ad
109 : else
110 : esum_ad = - dL_dm_ad + sources_ad + &
111 52348 : others_ad - d_turbulent_energy_dt_ad - dwork_dm_ad - dke_dt_ad - dpe_dt_ad - de_dt_ad
112 : end if
113 52348 : resid_ad = esum_ad ! convert back to auto_diff_real_star_order1
114 52348 : s% ergs_error(k) = -dm*dt*resid_ad%val ! save ergs_error before scaling
115 52348 : resid_ad = scal*resid_ad
116 52348 : residual = resid_ad%val
117 52348 : s% equ(i_dlnE_dt, k) = residual
118 :
119 : if (test_partials) then
120 : s% solver_test_partials_val = residual
121 : end if
122 52348 : call unpack_res18(s% species, resid_ad)
123 :
124 52348 : if (test_partials) then
125 : s% solver_test_partials_var = s% i_u
126 : s% solver_test_partials_dval_dx = d_d00(s% solver_test_partials_var)
127 : write(*,*) 'get1_energy_eqn', s% solver_test_partials_var
128 : if (eps_grav_form) write(*,*) 'eps_grav_form', eps_grav_form
129 : !if (.false. .and. s% solver_iter == s% solver_test_partials_iter_number) then
130 : if (.true.) then
131 : write(*,2) 'scal', k, scal
132 : write(*,2) 'residual', k, residual
133 : write(*,2) 'sources*scal', k, sources_ad%val*scal
134 : write(*,2) '-dL_dm*scal', k, -dL_dm_ad%val*scal
135 : write(*,2) '-d_turbulent_energy_dt*scal', k, -d_turbulent_energy_dt_ad%val*scal
136 : write(*,2) '-dwork_dm*scal', k, -dwork_dm_ad%val*scal
137 : write(*,2) '-dke_dt*scal', k, -dke_dt_ad%val*scal
138 : write(*,2) '-dpe_dt*scal', k, -dpe_dt_ad%val*scal
139 : write(*,2) 'gradT', k, s% gradT(k)
140 : write(*,2) 'opacity', k, s% opacity(k)
141 : write(*,2) 'logT', k, s% lnT(k)/ln10
142 : write(*,2) 'logRho', k, s% lnd(k)/ln10
143 : write(*,2) 'X', k, s% X(k)
144 : write(*,2) 'Z', k, s% Z(k)
145 : end if
146 : write(*,'(A)')
147 : end if
148 :
149 : contains
150 :
151 52348 : subroutine init
152 52348 : i_dlnE_dt = s% i_dlnE_dt
153 52348 : i_lum = s% i_lum
154 52348 : i_v = s% i_v
155 52348 : nz = s% nz
156 52348 : dt = s% dt
157 52348 : dm = s% dm(k)
158 : doing_op_split_burn = s% op_split_burn .and. &
159 52348 : s% T_start(k) >= s% op_split_burn_min_T
160 1936876 : d_dm1 = 0d0; d_d00 = 0d0; d_dp1 = 0d0
161 52348 : end subroutine init
162 :
163 52348 : subroutine setup_dwork_dm(ierr)
164 : integer, intent(out) :: ierr
165 : real(dp) :: dwork
166 : logical :: skip_P
167 : include 'formats'
168 52348 : ierr = 0
169 52348 : skip_P = eps_grav_form
170 52348 : if (s% using_velocity_time_centering .and. &
171 : s% use_P_d_1_div_rho_form_of_work_when_time_centering_velocity) then
172 : call eval_simple_PdV_work(s, k, skip_P, dwork_dm_ad, dwork, &
173 0 : d_dwork_dxa00, ierr)
174 0 : d_dwork_dxam1 = 0
175 0 : d_dwork_dxap1 = 0
176 0 : if (k == s% nz) then
177 0 : s% work_inward_at_center = pi4*pow2(s% r_center)*s% Peos_start(s% nz)*s% v_center
178 0 : if (is_bad(s% work_inward_at_center)) then
179 0 : write(*,2) 'work_inward_at_center', s% model_number, s% work_inward_at_center
180 0 : write(*,2) 'Peos_start', s% model_number, s% Peos_start(s% nz)
181 0 : write(*,2) 'v_center', s% model_number, s% v_center
182 0 : write(*,2) 'r_center', s% model_number, s% r_center
183 0 : call mesa_error(__FILE__,__LINE__,'setup_dwork_dm')
184 : end if
185 : end if
186 : else
187 : call eval_dwork(s, k, skip_P, dwork_dm_ad, dwork, &
188 52348 : d_dwork_dxam1, d_dwork_dxa00, d_dwork_dxap1, ierr)
189 : end if
190 52348 : if (ierr /= 0) then
191 0 : if (s% report_ierr) write(*,*) 'failed in setup_dwork_dm', k
192 0 : return
193 : end if
194 52348 : dwork_dm_ad = dwork_dm_ad/dm
195 : end subroutine setup_dwork_dm
196 :
197 52348 : subroutine setup_dL_dm(ierr)
198 : integer, intent(out) :: ierr
199 : type(auto_diff_real_star_order1) :: L00_ad, Lp1_ad
200 : real(dp) :: L_theta
201 : include 'formats'
202 52348 : ierr = 0
203 : if (s% using_velocity_time_centering .and. &
204 : s% include_L_in_velocity_time_centering &
205 52348 : .and. s% lnT(k)/ln10 <= s% max_logT_for_include_P_and_L_in_velocity_time_centering) then
206 0 : L_theta = s% L_theta_for_velocity_time_centering
207 : else
208 52348 : L_theta = 1d0
209 : end if
210 52348 : L00_ad = L_theta*wrap_L_00(s, k) + (1d0 - L_theta)*s% L_start(k)
211 52348 : Lp1_ad = wrap_L_p1(s, k)
212 52348 : if (k < s% nz) Lp1_ad = L_theta*Lp1_ad + (1d0 - L_theta)*s% L_start(k+1)
213 52348 : dL_dm_ad = (L00_ad - Lp1_ad)/dm
214 52348 : end subroutine setup_dL_dm
215 :
216 :
217 52348 : subroutine setup_sources_and_others(ierr) ! sources_ad, others_ad
218 : use hydro_rsp2, only: compute_Eq_cell
219 : use star_utils, only: get_face_weights
220 : use tdc_hydro, only: compute_tdc_Eq_div_w_face ! compute_Eq_cell
221 52348 : real(dp) :: alfa, beta
222 : integer, intent(out) :: ierr
223 : type(auto_diff_real_star_order1) :: &
224 : eps_nuc_ad, non_nuc_neu_ad, extra_heat_ad, Eq_ad, RTI_diffusion_ad, &
225 : v_00, v_p1, drag_force, drag_energy
226 : include 'formats'
227 52348 : ierr = 0
228 :
229 52348 : if (s% eps_nuc_factor == 0d0 .or. s% nonlocal_NiCo_decay_heat) then
230 0 : eps_nuc_ad = 0 ! get eps_nuc from extra_heat instead
231 52348 : else if (s% op_split_burn .and. s% T_start(k) >= s% op_split_burn_min_T) then
232 0 : eps_nuc_ad = 0d0
233 0 : eps_nuc_ad%val = s% burn_avg_epsnuc(k)
234 : else
235 52348 : eps_nuc_ad = 0d0
236 52348 : eps_nuc_ad%val = s% eps_nuc(k)
237 52348 : eps_nuc_ad%d1Array(i_lnd_00) = s% d_epsnuc_dlnd(k)
238 52348 : eps_nuc_ad%d1Array(i_lnT_00) = s% d_epsnuc_dlnT(k)
239 : end if
240 :
241 52348 : non_nuc_neu_ad = 0d0
242 : ! for reasons lost in the past, we always time center non_nuc_neu
243 : ! change that if you are feeling lucky.
244 52348 : non_nuc_neu_ad%val = 0.5d0*(s% non_nuc_neu_start(k) + s% non_nuc_neu(k))
245 52348 : non_nuc_neu_ad%d1Array(i_lnd_00) = 0.5d0*s% d_nonnucneu_dlnd(k)
246 52348 : non_nuc_neu_ad%d1Array(i_lnT_00) = 0.5d0*s% d_nonnucneu_dlnT(k)
247 :
248 52348 : extra_heat_ad = s% extra_heat(k)
249 :
250 : ! other = eps_WD_sedimentation + eps_diffusion + eps_pre_mix + eps_phase_separation
251 : ! no partials for any of these
252 52348 : others_ad = 0d0
253 52348 : if (s% do_element_diffusion) then
254 0 : if (s% do_WD_sedimentation_heating) then
255 0 : others_ad%val = others_ad%val + s% eps_WD_sedimentation(k)
256 0 : else if (s% do_diffusion_heating) then
257 0 : others_ad%val = others_ad%val + s% eps_diffusion(k)
258 : end if
259 : end if
260 52348 : if (s% do_conv_premix .and. s% do_premix_heating) &
261 0 : others_ad%val = others_ad%val + s% eps_pre_mix(k)
262 52348 : if (s% do_phase_separation .and. s% do_phase_separation_heating) &
263 0 : others_ad%val = others_ad%val + s% eps_phase_separation(k)
264 :
265 52348 : Eq_ad = 0d0
266 52348 : if (s% RSP2_flag) then
267 0 : Eq_ad = s% Eq_ad(k) ! compute_Eq_cell(s, k, ierr)
268 0 : if (ierr /= 0) return
269 52348 : else if (s% alpha_TDC_DampM >0d0 .and. s% MLT_option == 'TDC' .and. &
270 : s% TDC_include_eturb_in_energy_equation) then ! not checking for v or u flag.
271 : !Eq_ad = compute_tdc_Eq_cell(s, k, ierr) ! safe to just recompute
272 0 : if (k == 1) then
273 0 : Eq_ad = compute_tdc_Eq_div_w_face(s, k, ierr)*(s% mlt_vc_ad(k)/sqrt_2_div_3)
274 : else
275 0 : call get_face_weights(s, k, alfa, beta)
276 : Eq_ad = alfa*compute_tdc_Eq_div_w_face(s, k, ierr)*(s% mlt_vc_ad(k)/sqrt_2_div_3) + &
277 0 : beta*shift_m1(compute_tdc_Eq_div_w_face(s, k-1, ierr))*(shift_m1(s% mlt_vc_ad(k-1))/sqrt_2_div_3)
278 : end if
279 0 : if (ierr /= 0) return
280 : end if
281 :
282 52348 : call setup_RTI_diffusion(RTI_diffusion_ad)
283 :
284 52348 : drag_energy = 0d0
285 52348 : s% FdotV_drag_energy(k) = 0
286 52348 : if (k /= s% nz) then
287 : if ((s% q(k) > s% min_q_for_drag) .and. &
288 52304 : (s% drag_coefficient > 0) .and. &
289 : s% use_drag_energy) then
290 0 : v_00 = wrap_v_00(s,k)
291 0 : drag_force = s% drag_coefficient*v_00/s% dt
292 0 : drag_energy = 0.5d0*v_00*drag_force
293 0 : s% FdotV_drag_energy(k) = drag_energy%val
294 : ! drag energy for outer half-cell. the 0.5d0 is for dm/2
295 : end if
296 : if ((s% q(k+1) > s% min_q_for_drag) .and. &
297 52304 : (s% drag_coefficient > 0) .and. &
298 : s% use_drag_energy) then
299 0 : v_p1 = wrap_v_p1(s,k)
300 0 : drag_force = s% drag_coefficient*v_p1/s% dt
301 0 : drag_energy = drag_energy + 0.5d0*v_p1*drag_force
302 0 : s% FdotV_drag_energy(k) = drag_energy%val
303 : ! drag energy for inner half-cell. the 0.5d0 is for dm/2
304 : end if
305 : end if
306 :
307 52348 : sources_ad = eps_nuc_ad - non_nuc_neu_ad + extra_heat_ad + Eq_ad + RTI_diffusion_ad + drag_energy
308 :
309 52348 : sources_ad%val = sources_ad%val + s% irradiation_heat(k)
310 :
311 52348 : if (s% mstar_dot /= 0d0) sources_ad%val = sources_ad%val + s% eps_mdot(k)
312 :
313 52348 : end subroutine setup_sources_and_others
314 :
315 52348 : subroutine setup_RTI_diffusion(diffusion_eps_ad)
316 : type(auto_diff_real_star_order1), intent(out) :: diffusion_eps_ad
317 52348 : real(dp) :: diffusion_factor, emin_start, sigp1, sig00
318 : logical :: do_diffusion
319 : type(auto_diff_real_star_order1) :: &
320 : e_m1, e_00, e_p1, diffusion_eps_in, diffusion_eps_out
321 : include 'formats'
322 52348 : diffusion_factor = s% dedt_RTI_diffusion_factor
323 52348 : do_diffusion = s% RTI_flag .and. diffusion_factor > 0d0
324 : if (.not. do_diffusion) then
325 52348 : diffusion_eps_ad = 0d0
326 : else
327 0 : if (k < s% nz) then
328 0 : if (s% alpha_RTI(k) > 1d-10 .and. k > 1) then
329 : emin_start = min( &
330 0 : s% energy_start(k+1), s% energy_start(k), s% energy_start(k-1))
331 0 : if (emin_start < 5d0*s% RTI_energy_floor) then
332 : diffusion_factor = diffusion_factor* &
333 0 : (1d0 + (5d0*s% RTI_energy_floor - emin_start)/emin_start)
334 : end if
335 : end if
336 0 : sigp1 = diffusion_factor*s% sig_RTI(k+1)
337 0 : e_p1 = wrap_e_p1(s,k)
338 : else
339 0 : sigp1 = 0
340 0 : e_p1 = 0d0
341 : end if
342 0 : if (k > 1) then
343 0 : sig00 = diffusion_factor*s% sig_RTI(k)
344 0 : e_m1 = wrap_e_m1(s,k)
345 : else
346 0 : sig00 = 0
347 0 : e_m1 = 0
348 : end if
349 0 : e_00 = wrap_e_00(s,k)
350 0 : diffusion_eps_in = sigp1*(e_p1 - e_00)/dm
351 0 : diffusion_eps_out = sig00*(e_00 - e_m1)/dm
352 0 : diffusion_eps_ad = diffusion_eps_in - diffusion_eps_out
353 : end if
354 52348 : s% dedt_RTI(k) = diffusion_eps_ad%val
355 52348 : end subroutine setup_RTI_diffusion
356 :
357 52348 : subroutine setup_d_turbulent_energy_dt(ierr)
358 : use star_utils, only: get_face_weights
359 : use const_def, only: sqrt_2_div_3
360 : integer, intent(out) :: ierr
361 : type(auto_diff_real_star_order1) :: TDC_eturb_cell
362 52348 : real (dp) :: TDC_eturb_cell_start, alfa, beta
363 : include 'formats'
364 52348 : ierr = 0
365 52348 : if (s% RSP2_flag) then
366 0 : d_turbulent_energy_dt_ad = (wrap_etrb_00(s,k) - get_etrb_start(s,k))/dt
367 52348 : else if (s% mlt_vc_old(k) > 0d0 .and. s% MLT_option == 'TDC' .and. &
368 : s% TDC_include_eturb_in_energy_equation) then
369 : ! write a wrapper for this.
370 0 : if (k == 1) then
371 0 : TDC_eturb_cell_start = pow2(s% mlt_vc_old(k)/sqrt_2_div_3)
372 0 : TDC_eturb_cell = pow2(s% mlt_vc(k)/sqrt_2_div_3)
373 : else
374 0 : call get_face_weights(s, k, alfa, beta)
375 : TDC_eturb_cell_start = alfa*pow2(s% mlt_vc_old(k)/sqrt_2_div_3) + &
376 0 : beta*pow2(s% mlt_vc_old(k-1)/sqrt_2_div_3)
377 : TDC_eturb_cell = alfa*pow2(s% mlt_vc_ad(k)/sqrt_2_div_3) + &
378 0 : beta*pow2(shift_m1(s% mlt_vc_ad(k-1))/sqrt_2_div_3)
379 : end if
380 0 : d_turbulent_energy_dt_ad = (TDC_eturb_cell - TDC_eturb_cell_start) / dt
381 : else
382 52348 : d_turbulent_energy_dt_ad = 0d0
383 : end if
384 52348 : s% detrbdt(k) = d_turbulent_energy_dt_ad%val
385 52348 : end subroutine setup_d_turbulent_energy_dt
386 :
387 52348 : subroutine setup_eps_grav(ierr)
388 : integer, intent(out) :: ierr
389 : include 'formats'
390 52348 : ierr = 0
391 :
392 52348 : if (s% u_flag) then ! for now, assume u_flag means no eps_grav
393 0 : eps_grav_form = .false.
394 0 : return
395 : end if
396 :
397 : ! value from checking s% energy_eqn_option in hydro_eqns.f90
398 52348 : eps_grav_form = s% eps_grav_form_for_energy_eqn
399 :
400 52348 : if (.not. eps_grav_form) then ! check if want it true
401 52348 : if (s% doing_relax .and. s% no_dedt_form_during_relax) eps_grav_form = .true.
402 : end if
403 :
404 52348 : if (eps_grav_form) then
405 0 : if (s% RSP2_flag) then
406 0 : call mesa_error(__FILE__,__LINE__,'cannot use eps_grav with et yet. fix energy eqn.')
407 : end if
408 0 : call eval_eps_grav_and_partials(s, k, ierr) ! get eps_grav info
409 0 : if (ierr /= 0) then
410 0 : if (s% report_ierr) write(*,2) 'failed in eval_eps_grav_and_partials', k
411 0 : return
412 : end if
413 0 : eps_grav_ad = s% eps_grav_ad(k)
414 : end if
415 :
416 52348 : end subroutine setup_eps_grav
417 :
418 52348 : subroutine setup_de_dt_and_friends(ierr)
419 : use star_utils, only: get_dke_dt_dpe_dt
420 : integer, intent(out) :: ierr
421 : real(dp) :: dke_dt, d_dkedt_dv00, d_dkedt_dvp1, &
422 : dpe_dt, d_dpedt_dlnR00, d_dpedt_dlnRp1, &
423 52348 : de_dt, d_de_dt_dlnd, d_de_dt_dlnT
424 : include 'formats'
425 52348 : ierr = 0
426 :
427 52348 : dke_dt = 0d0; d_dkedt_dv00 = 0d0; d_dkedt_dvp1 = 0d0
428 52348 : dpe_dt = 0d0; d_dpedt_dlnR00 = 0d0; d_dpedt_dlnRp1 = 0d0
429 52348 : de_dt = 0d0; d_de_dt_dlnd = 0d0; d_de_dt_dlnT = 0d0
430 :
431 52348 : if (.not. eps_grav_form) then
432 :
433 52348 : de_dt = (s% energy(k) - s% energy_start(k))/dt
434 52348 : d_de_dt_dlnd = s% dE_dRho_for_partials(k)*s% rho(k)/dt
435 52348 : d_de_dt_dlnT = s% Cv_for_partials(k)*s% T(k)/dt
436 52348 : de_dt_ad = 0d0
437 52348 : de_dt_ad%val = de_dt
438 52348 : de_dt_ad%d1Array(i_lnd_00) = d_de_dt_dlnd
439 52348 : de_dt_ad%d1Array(i_lnT_00) = d_de_dt_dlnT
440 :
441 : call get_dke_dt_dpe_dt(s, k, dt, &
442 : dke_dt, d_dkedt_dv00, d_dkedt_dvp1, &
443 52348 : dpe_dt, d_dpedt_dlnR00, d_dpedt_dlnRp1, ierr)
444 52348 : if (ierr /= 0) then
445 0 : if (s% report_ierr) write(*,2) 'failed in get_dke_dt_dpe_dt', k
446 : return
447 : end if
448 52348 : dke_dt_ad = 0d0
449 52348 : dke_dt_ad%val = dke_dt
450 52348 : dke_dt_ad%d1Array(i_v_00) = d_dkedt_dv00
451 52348 : dke_dt_ad%d1Array(i_v_p1) = d_dkedt_dvp1
452 :
453 52348 : dpe_dt_ad = 0d0
454 52348 : dpe_dt_ad%val = dpe_dt
455 52348 : dpe_dt_ad%d1Array(i_lnR_00) = d_dpedt_dlnR00
456 52348 : dpe_dt_ad%d1Array(i_lnR_p1) = d_dpedt_dlnRp1
457 :
458 : end if
459 :
460 52348 : s% dkedt(k) = dke_dt
461 52348 : s% dpedt(k) = dpe_dt
462 52348 : s% dkedt(k) = dke_dt
463 52348 : s% dedt(k) = de_dt
464 :
465 52348 : end subroutine setup_de_dt_and_friends
466 :
467 52348 : subroutine unpack_res18(species,res18)
468 52348 : use star_utils, only: save_eqn_dxa_partials, unpack_residual_partials
469 : type(auto_diff_real_star_order1) :: res18
470 : integer, intent(in) :: species
471 52348 : real(dp) :: dequ
472 : integer :: j
473 1413396 : real(dp), dimension(species) :: dxam1, dxa00, dxap1
474 : logical, parameter :: checking = .true.
475 : include 'formats'
476 :
477 : ! do partials wrt composition
478 1413396 : dxam1 = 0d0; dxa00 = 0d0; dxap1 = 0d0
479 52348 : if (.not. (s% nonlocal_NiCo_decay_heat .or. doing_op_split_burn)) then
480 52348 : if (do_chem .and. s% dxdt_nuc_factor > 0d0) then
481 471132 : do j=1,s% species
482 418784 : dequ = scal*s% d_epsnuc_dx(j,k)
483 418784 : if (checking) call check_dequ(dequ,'d_epsnuc_dx')
484 471132 : dxa00(j) = dxa00(j) + dequ
485 : end do
486 : end if
487 : end if
488 :
489 52348 : if (.not. eps_grav_form) then
490 471132 : do j=1,s% species
491 418784 : dequ = -scal*(s%energy(k)/dt)*s% dlnE_dxa_for_partials(j,k)
492 418784 : if (checking) call check_dequ(dequ,'dlnE_dxa_for_partials')
493 471132 : dxa00(j) = dxa00(j) + dequ
494 : end do
495 0 : else if (do_chem .and. (.not. doing_op_split_burn) .and. &
496 : (s% dxdt_nuc_factor > 0d0 .or. s% mix_factor > 0d0)) then
497 0 : do j=1,s% species
498 0 : dequ = scal*s% d_eps_grav_dx(j,k)
499 0 : if (checking) call check_dequ(dequ,'d_eps_grav_dx')
500 0 : dxa00(j) = dxa00(j) + dequ
501 : end do
502 : end if
503 :
504 471132 : do j=1,s% species
505 418784 : dequ = -scal*d_dwork_dxa00(j)/dm
506 418784 : if (checking) call check_dequ(dequ,'d_dwork_dxa00')
507 471132 : dxa00(j) = dxa00(j) + dequ
508 : end do
509 52348 : if (k > 1) then
510 470736 : do j=1,s% species
511 418432 : dequ = -scal*d_dwork_dxam1(j)/dm
512 418432 : if (checking) call check_dequ(dequ,'d_dwork_dxam1')
513 470736 : dxam1(j) = dxam1(j) + dequ
514 : end do
515 : end if
516 52348 : if (k < nz) then
517 470736 : do j=1,s% species
518 418432 : dequ = -scal*d_dwork_dxap1(j)/dm
519 418432 : if (checking) call check_dequ(dequ,'d_dwork_dxap1')
520 470736 : dxap1(j) = dxap1(j) + dequ
521 : end do
522 : end if
523 :
524 : call save_eqn_dxa_partials(&
525 52348 : s, k, nvar, i_dlnE_dt, species, dxam1, dxa00, dxap1, 'get1_energy_eqn', ierr)
526 :
527 : call unpack_residual_partials(s, k, nvar, i_dlnE_dt, &
528 52348 : res18, d_dm1, d_d00, d_dp1)
529 :
530 52348 : end subroutine unpack_res18
531 :
532 2093216 : subroutine check_dequ(dequ, str)
533 : real(dp), intent(in) :: dequ
534 : character (len=*), intent(in) :: str
535 : include 'formats'
536 2093216 : if (is_bad(dequ)) then
537 0 : !$omp critical (hydro_energy_crit2)
538 0 : ierr = -1
539 0 : if (s% report_ierr) then
540 0 : write(*,2) 'get1_energy_eqn: bad ' // trim(str), k, dequ
541 : end if
542 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'get1_energy_eqn')
543 : !$omp end critical (hydro_energy_crit2)
544 0 : return
545 : end if
546 52348 : end subroutine check_dequ
547 :
548 : subroutine unpack1(j, dvar_m1, dvar_00, dvar_p1)
549 : integer, intent(in) :: j
550 : real(dp), intent(in) :: dvar_m1, dvar_00, dvar_p1
551 : d_dm1(j) = dvar_m1
552 : d_d00(j) = dvar_00
553 : d_dp1(j) = dvar_p1
554 : end subroutine unpack1
555 :
556 : end subroutine get1_energy_eqn
557 :
558 :
559 52348 : subroutine eval_dwork(s, k, skip_P, dwork_ad, dwork, &
560 52348 : d_dwork_dxam1, d_dwork_dxa00, d_dwork_dxap1, ierr)
561 : use accurate_sum_auto_diff_star_order1
562 : use auto_diff_support
563 : use star_utils, only: calc_Ptot_ad_tw
564 : type (star_info), pointer :: s
565 : integer, intent(in) :: k
566 : logical, intent(in) :: skip_P
567 : type(auto_diff_real_star_order1), intent(out) :: dwork_ad
568 : real(dp), intent(out) :: dwork
569 : real(dp), intent(out), dimension(s% species) :: &
570 : d_dwork_dxam1, d_dwork_dxa00, d_dwork_dxap1
571 : integer, intent(out) :: ierr
572 :
573 52348 : real(dp) :: work_00, work_p1
574 : real(dp), dimension(s% species) :: &
575 889916 : d_work_00_dxa00, d_work_00_dxam1, &
576 889916 : d_work_p1_dxap1, d_work_p1_dxa00
577 : type(auto_diff_real_star_order1) :: work_00_ad, work_p1_ad
578 : logical :: test_partials
579 : integer :: j
580 : include 'formats'
581 : ierr = 0
582 :
583 : call eval1_work(s, k, skip_P, &
584 52348 : work_00_ad, work_00, d_work_00_dxa00, d_work_00_dxam1, ierr)
585 52348 : if (ierr /= 0) return
586 : call eval1_work(s, k+1, skip_P, &
587 52348 : work_p1_ad, work_p1, d_work_p1_dxap1, d_work_p1_dxa00, ierr)
588 52348 : if (ierr /= 0) return
589 52348 : work_p1_ad = shift_p1(work_p1_ad) ! shift the partials
590 52348 : dwork_ad = work_00_ad - work_p1_ad
591 52348 : dwork = dwork_ad%val
592 471132 : do j=1,s% species
593 418784 : d_dwork_dxam1(j) = d_work_00_dxam1(j)
594 418784 : d_dwork_dxa00(j) = d_work_00_dxa00(j) - d_work_p1_dxa00(j)
595 471132 : d_dwork_dxap1(j) = -d_work_p1_dxap1(j)
596 : end do
597 :
598 : !test_partials = (k == s% solver_test_partials_k)
599 52348 : test_partials = .false.
600 :
601 : if (test_partials) then
602 : s% solver_test_partials_val = 0
603 : s% solver_test_partials_var = 0
604 : s% solver_test_partials_dval_dx = 0
605 : write(*,*) 'eval_dwork', s% solver_test_partials_var
606 : end if
607 :
608 52348 : end subroutine eval_dwork
609 :
610 :
611 : ! ergs/s at face(k)
612 104696 : subroutine eval1_work(s, k, skip_Peos, &
613 104696 : work_ad, work, d_work_dxa00, d_work_dxam1, ierr)
614 52348 : use star_utils, only: get_Pvsc_ad, calc_Ptrb_ad_tw, get_rho_face
615 : use accurate_sum_auto_diff_star_order1
616 : use auto_diff_support
617 : type (star_info), pointer :: s
618 : integer, intent(in) :: k
619 : logical, intent(in) :: skip_Peos
620 : type(auto_diff_real_star_order1), intent(out) :: work_ad
621 : real(dp), intent(out) :: work
622 : real(dp), dimension(s% species), intent(out) :: &
623 : d_work_dxa00, d_work_dxam1
624 : integer, intent(out) :: ierr
625 209392 : real(dp) :: alfa, beta, P_theta, Av_face
626 1779832 : real(dp), dimension(s% species) :: d_Pface_dxa00, d_Pface_dxam1
627 : type(auto_diff_real_star_order1) :: &
628 : P_face_ad, A_times_v_face_ad, mlt_Pturb_ad, &
629 : PtrbR_ad, PtrbL_ad, PvscL_ad, PvscR_ad, Ptrb_div_etrb, PL_ad, PR_ad, &
630 : Peos_ad, Ptrb_ad, Pvsc_ad, extra_P
631 : logical :: test_partials
632 : integer :: j
633 : include 'formats'
634 104696 : ierr = 0
635 :
636 942264 : d_work_dxa00 = 0d0
637 942264 : d_work_dxam1 = 0d0
638 104696 : if (k > s% nz .or. (s% dt <= 0d0 .and. .not. (s% v_flag .or. s% u_flag))) then
639 44 : work_ad = 0d0
640 44 : if (k == s% nz+1) then
641 44 : work = pi4*pow2(s% r_center)*s% Peos_start(s% nz)*s% v_center
642 44 : s% work_inward_at_center = work
643 44 : if (is_bad(work)) then
644 0 : write(*,2) 'work_inward_at_center', s% model_number, work
645 0 : write(*,2) 'Peos_start', s% model_number, s% Peos_start(s% nz)
646 0 : write(*,2) 'v_center', s% model_number, s% v_center
647 0 : write(*,2) 'r_center', s% model_number, s% r_center
648 0 : call mesa_error(__FILE__,__LINE__,'eval1_work')
649 : end if
650 : end if
651 44 : work_ad%val = work
652 44 : return
653 : end if
654 :
655 104652 : call eval1_A_times_v_face_ad(s, k, A_times_v_face_ad, ierr)
656 104652 : if (ierr /= 0) return
657 :
658 104652 : if (k > 1) then
659 104608 : alfa = s% dq(k-1)/(s% dq(k-1) + s% dq(k))
660 : else
661 44 : alfa = 1d0
662 : end if
663 104652 : beta = 1d0 - alfa
664 :
665 : if (s% using_velocity_time_centering .and. &
666 104652 : s% include_P_in_velocity_time_centering .and. &
667 : s% lnT(k)/ln10 <= s% max_logT_for_include_P_and_L_in_velocity_time_centering) then
668 0 : P_theta = s% P_theta_for_velocity_time_centering
669 : else
670 104652 : P_theta = 1d0 ! try 1 - q(k)
671 : end if
672 :
673 104652 : if (s% u_flag) then
674 0 : P_face_ad = P_theta*s% P_face_ad(k) + (1d0-P_theta)*s% P_face_start(k)
675 0 : d_Pface_dxa00 = 0d0
676 0 : d_Pface_dxam1 = 0d0
677 : else ! set P_ad
678 941868 : d_Pface_dxa00 = 0d0
679 941868 : d_Pface_dxam1 = 0d0
680 104652 : if (skip_Peos) then
681 0 : Peos_ad = 0d0
682 : else
683 104652 : if (k > 1) then
684 104608 : PR_ad = P_theta*wrap_Peos_m1(s,k) + (1d0-P_theta)*s% Peos_start(k-1)
685 : else
686 44 : PR_ad = 0d0
687 : end if
688 104652 : PL_ad = P_theta*wrap_Peos_00(s,k) + (1d0-P_theta)*s% Peos_start(k)
689 104652 : Peos_ad = alfa*PL_ad + beta*PR_ad
690 104652 : if (k > 1) then
691 941472 : do j=1,s% species
692 : d_Pface_dxa00(j) = &
693 941472 : alfa*s% dlnPeos_dxa_for_partials(j,k)*P_theta*s% Peos(k)
694 : end do
695 941472 : do j=1,s% species
696 : d_Pface_dxam1(j) = &
697 941472 : beta*s% dlnPeos_dxa_for_partials(j,k-1)*P_theta*s% Peos(k-1)
698 : end do
699 : else ! k == 1
700 396 : do j=1,s% species
701 : d_Pface_dxa00(j) = &
702 396 : s% dlnPeos_dxa_for_partials(j,k)*P_theta*s% Peos(k)
703 : end do
704 : end if
705 : end if
706 :
707 : ! set Pvsc_ad
708 104652 : if (.not. s% use_Pvsc_art_visc) then
709 104652 : Pvsc_ad = 0d0
710 : else
711 0 : if (k > 1) then
712 0 : call get_Pvsc_ad(s, k-1, PvscR_ad, ierr)
713 0 : if (ierr /= 0) return
714 0 : PvscR_ad = shift_m1(PvscR_ad)
715 0 : if (s% include_P_in_velocity_time_centering .and. &
716 : s% lnT(k)/ln10 <= s% max_logT_for_include_P_and_L_in_velocity_time_centering) &
717 0 : PvscR_ad = 0.5d0*(PvscR_ad + s% Pvsc_start(k-1))
718 : else
719 0 : PvscR_ad = 0d0
720 : end if
721 0 : call get_Pvsc_ad(s, k, PvscL_ad, ierr)
722 0 : if (ierr /= 0) return
723 0 : if (s% include_P_in_velocity_time_centering .and. &
724 : s% lnT(k)/ln10 <= s% max_logT_for_include_P_and_L_in_velocity_time_centering) &
725 0 : PvscL_ad = 0.5d0*(PvscL_ad + s% Pvsc_start(k))
726 0 : Pvsc_ad = alfa*PvscL_ad + beta*PvscR_ad
727 : end if
728 :
729 : ! set Ptrb_ad
730 104652 : if (.not. s% RSP2_flag) then
731 104652 : Ptrb_ad = 0d0
732 : else
733 0 : if (k > 1) then
734 0 : call calc_Ptrb_ad_tw(s, k-1, PtrbR_ad, Ptrb_div_etrb, ierr)
735 0 : if (ierr /= 0) return
736 0 : PtrbR_ad = shift_m1(PtrbR_ad)
737 : else
738 0 : PtrbR_ad = 0d0
739 : end if
740 0 : call calc_Ptrb_ad_tw(s, k, PtrbL_ad, Ptrb_div_etrb, ierr)
741 0 : if (ierr /= 0) return
742 0 : Ptrb_ad = alfa*PtrbL_ad + beta*PtrbR_ad
743 : end if
744 :
745 : ! set extra_P
746 104652 : if (.not. s% use_other_pressure) then
747 104652 : extra_P = 0d0
748 0 : else if (k > 1) then
749 : ! my_val_m1 = shift_m1(get_my_val(s,k-1)) for use in terms going into equation at k
750 0 : extra_P = alfa*s% extra_pressure(k) + beta * shift_m1(s%extra_pressure(k-1))
751 : else
752 0 : extra_P = s% extra_pressure(k)
753 : end if
754 :
755 : ! set mlt_Pturb_ad
756 104652 : mlt_Pturb_ad = 0d0
757 104652 : if (s% mlt_Pturb_factor > 0d0 .and. s% mlt_vc_old(k) > 0d0) &
758 0 : mlt_Pturb_ad = s% mlt_Pturb_factor*pow2(s% mlt_vc_old(k))*get_rho_face(s,k)/3d0
759 :
760 104652 : P_face_ad = Peos_ad + Pvsc_ad + Ptrb_ad + mlt_Pturb_ad + extra_P
761 :
762 : end if
763 :
764 104652 : work_ad = A_times_v_face_ad*P_face_ad
765 104652 : work = work_ad%val
766 :
767 104652 : if (k == 1) s% work_outward_at_surface = work
768 :
769 104652 : Av_face = A_times_v_face_ad%val
770 941868 : do j=1,s% species
771 837216 : d_work_dxa00(j) = Av_face*d_Pface_dxa00(j)
772 941868 : d_work_dxam1(j) = Av_face*d_Pface_dxam1(j)
773 : end do
774 :
775 : !test_partials = (k == s% solver_test_partials_k)
776 104652 : test_partials = .false.
777 :
778 : if (test_partials) then
779 : s% solver_test_partials_val = 0
780 : s% solver_test_partials_var = 0
781 : s% solver_test_partials_dval_dx = 0
782 : write(*,*) 'eval1_work', s% solver_test_partials_var
783 : end if
784 :
785 104696 : end subroutine eval1_work
786 :
787 :
788 104652 : subroutine eval1_A_times_v_face_ad(s, k, A_times_v_face_ad, ierr)
789 104696 : use star_utils, only: get_area_info_opt_time_center
790 : type (star_info), pointer :: s
791 : integer, intent(in) :: k
792 : type(auto_diff_real_star_order1), intent(out) :: A_times_v_face_ad
793 : integer, intent(out) :: ierr
794 : type(auto_diff_real_star_order1) :: A_ad, inv_R2, u_face_ad
795 : include 'formats'
796 :
797 : ierr = 0
798 104652 : call get_area_info_opt_time_center(s, k, A_ad, inv_R2, ierr)
799 104652 : if (ierr /= 0) return
800 :
801 104652 : u_face_ad = 0d0
802 104652 : if (s% v_flag) then
803 0 : u_face_ad%val = s% vc(k)
804 0 : u_face_ad%d1Array(i_v_00) = s% d_vc_dv
805 104652 : else if (s% u_flag) then
806 0 : u_face_ad = s% u_face_ad(k)
807 0 : if (s% using_velocity_time_centering) &
808 0 : u_face_ad = 0.5d0*(u_face_ad + s% u_face_start(k))
809 104652 : else if (s% using_velocity_time_centering) then
810 0 : u_face_ad%val = 0.5d0*(s% r(k) - s% r_start(k))/s% dt
811 0 : u_face_ad%d1Array(i_lnR_00) = 0.5d0*s% r(k)/s% dt
812 : else
813 104652 : u_face_ad%val = (s% r(k) - s% r_start(k))/s% dt
814 104652 : u_face_ad%d1Array(i_lnR_00) = s% r(k)/s% dt
815 : end if
816 :
817 104652 : A_times_v_face_ad = A_ad*u_face_ad
818 :
819 104652 : end subroutine eval1_A_times_v_face_ad
820 :
821 :
822 0 : subroutine eval_simple_PdV_work( &
823 0 : s, k, skip_P, dwork_ad, dwork, d_dwork_dxa00, ierr)
824 104652 : use accurate_sum_auto_diff_star_order1
825 : use auto_diff_support
826 : use star_utils, only: calc_Ptot_ad_tw
827 : type (star_info), pointer :: s
828 : integer, intent(in) :: k
829 : logical, intent(in) :: skip_P
830 : type(auto_diff_real_star_order1), intent(out) :: dwork_ad
831 : real(dp), intent(out) :: dwork
832 : real(dp), intent(out), dimension(s% species) :: d_dwork_dxa00
833 : integer, intent(out) :: ierr
834 :
835 : type(auto_diff_real_star_order1) :: &
836 : Av_face00_ad, Av_facep1_ad, Ptot_ad, dV
837 0 : real(dp), dimension(s% species) :: d_Ptot_dxa
838 0 : real(dp) :: Av_face00, Av_facep1
839 : logical :: include_mlt_Pturb
840 : integer :: j
841 :
842 : include 'formats'
843 : ierr = 0
844 :
845 : ! dV = 1/rho - 1/rho_start
846 0 : call eval1_A_times_v_face_ad(s, k, Av_face00_ad, ierr)
847 0 : if (ierr /= 0) return
848 0 : if (k < s% nz) then
849 0 : call eval1_A_times_v_face_ad(s, k+1, Av_facep1_ad, ierr)
850 0 : if (ierr /= 0) return
851 0 : Av_facep1_ad = shift_p1(Av_facep1_ad)
852 : else
853 0 : Av_facep1_ad = 0d0
854 0 : Av_facep1_ad%val = 4*pi*pow2(s% r_center)*s% v_center
855 : end if
856 0 : Av_face00 = Av_face00_ad%val
857 0 : Av_facep1 = Av_facep1_ad%val
858 0 : dV = Av_face00_ad - Av_facep1_ad
859 :
860 : include_mlt_Pturb = s% mlt_Pturb_factor > 0d0 &
861 0 : .and. s% mlt_vc_old(k) > 0d0 .and. k > 1
862 :
863 : call calc_Ptot_ad_tw( &
864 0 : s, k, skip_P, .not. include_mlt_Pturb, Ptot_ad, d_Ptot_dxa, ierr)
865 0 : if (ierr /= 0) return
866 :
867 0 : do j=1,s% species
868 0 : d_dwork_dxa00(j) = d_Ptot_dxa(j)*(Av_face00 - Av_facep1)
869 : end do
870 0 : if (k == 1) s% work_outward_at_surface = Ptot_ad%val*Av_face00
871 :
872 0 : dwork_ad = Ptot_ad*dV
873 0 : dwork = dwork_ad%val
874 :
875 0 : end subroutine eval_simple_PdV_work
876 :
877 : end module hydro_energy
|
