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