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