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