Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2012-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_eqns
21 :
22 : use star_private_def
23 : use const_def, only: dp, pi, ln10, two_thirds
24 : use star_utils, only: em1, e00, ep1
25 : use utils_lib, only: mesa_error, is_bad
26 : use auto_diff
27 : use auto_diff_support
28 :
29 : implicit none
30 :
31 : private
32 : public :: eval_equ
33 :
34 : contains
35 :
36 44 : subroutine eval_equ(s, nvar, ierr)
37 : type (star_info), pointer :: s
38 : integer, intent(in) :: nvar
39 : integer, intent(out) :: ierr
40 44 : call eval_equ_for_solver(s, nvar, 1, s% nz, ierr)
41 44 : end subroutine eval_equ
42 :
43 :
44 44 : subroutine eval_equ_for_solver(s, nvar, nzlo, nzhi, ierr)
45 : use chem_def
46 : use utils_lib, only: set_nan
47 : use mesh_functions
48 : use hydro_riemann, only: do_uface_and_Pface, do1_Riemann_momentum_eqn
49 : use hydro_momentum, only: do1_momentum_eqn, do1_radius_eqn
50 : use hydro_chem_eqns, only: do_chem_eqns, do1_chem_eqns
51 : use hydro_energy, only: do1_energy_eqn
52 : use hydro_temperature, only: do1_dlnT_dm_eqn
53 : use hydro_rsp2, only: do1_turbulent_energy_eqn, do1_rsp2_L_eqn, do1_rsp2_Hp_eqn
54 : use hydro_alpha_rti_eqns, only: do1_dalpha_RTI_dt_eqn
55 : use eps_grav, only: zero_eps_grav_and_partials
56 : use profile, only: do_save_profiles
57 : use star_utils, only: show_matrix, &
58 : no_extra_profile_columns, no_data_for_extra_profile_columns
59 :
60 : type (star_info), pointer :: s
61 : integer, intent(in) :: nvar, nzlo, nzhi
62 : integer, intent(out) :: ierr
63 :
64 : integer :: &
65 : i_dv_dt, i_du_dt, i_equL, i_dlnd_dt, i_dlnE_dt, i_dlnR_dt, &
66 : i_dalpha_RTI_dt, i_equ_w_div_wc, i_dj_rot_dt, i_detrb_dt, &
67 : i, k, j, nvar_hydro, nz, op_err
68 : integer :: &
69 : i_lnd, i_lnR, i_lnT, i_lum, i_v, i_u, i_w_div_wc, i_j_rot, &
70 : i_alpha_RTI, i_xh1, i_xhe4, species
71 44 : real(dp) :: L_phot_old
72 : real(dp), dimension(:), pointer :: &
73 44 : L, lnR, lnP, lnT, energy
74 : logical :: v_flag, u_flag, dump_for_debug, &
75 : do_chem, do_mix, do_dlnd_dt, do_dv_dt, do_du_dt, do_dlnR_dt, &
76 : do_alpha_RTI, do_w_div_wc, do_j_rot, do_dlnE_dt, do_equL, do_detrb_dt
77 :
78 : include 'formats'
79 :
80 44 : ierr = 0
81 :
82 44 : if (s% u_flag .and. s% use_mass_corrections) &
83 0 : call mesa_error(__FILE__,__LINE__,'use_mass_corrections dP not supported with u_flag true')
84 :
85 44 : if (s% u_flag) then
86 0 : call do_uface_and_Pface(s,ierr)
87 0 : if (ierr /= 0) then
88 0 : if (len_trim(s% retry_message) == 0) s% retry_message = 'do_uface_and_Pface failed'
89 0 : if (s% report_ierr) write(*,*) 'ierr from do_uface_and_Pface'
90 0 : return
91 : end if
92 : end if
93 :
94 44 : dump_for_debug = .false.
95 : !dump_for_debug = .true.
96 :
97 44 : do_mix = s% do_mix
98 44 : do_chem = (do_mix .or. s% do_burn)
99 :
100 44 : call unpack
101 :
102 : ! set flags indicating the variables currently in use
103 44 : do_dlnd_dt = (i_dlnd_dt > 0 .and. i_dlnd_dt <= nvar)
104 44 : do_dv_dt = (i_dv_dt > 0 .and. i_dv_dt <= nvar)
105 44 : do_du_dt = (i_du_dt > 0 .and. i_du_dt <= nvar)
106 44 : do_dlnR_dt = (i_dlnR_dt > 0 .and. i_dlnR_dt <= nvar)
107 44 : do_alpha_RTI = (i_alpha_RTI > 0 .and. i_alpha_RTI <= nvar)
108 44 : do_w_div_wc = (i_w_div_wc > 0 .and. i_w_div_wc <= nvar)
109 44 : do_j_rot = (i_j_rot > 0 .and. i_j_rot <= nvar)
110 44 : do_dlnE_dt = (i_dlnE_dt > 0 .and. i_dlnE_dt <= nvar)
111 44 : do_equL = (i_equL > 0 .and. i_equL <= nvar)
112 44 : do_detrb_dt = (i_detrb_dt > 0 .and. i_detrb_dt <= nvar)
113 :
114 44 : if (s% fill_arrays_with_NaNs) call set_nan(s% equ1)
115 :
116 : ! select form of energy equation
117 44 : if (trim(s% energy_eqn_option) == 'eps_grav') then
118 0 : s% eps_grav_form_for_energy_eqn = .true.
119 44 : else if (trim(s% energy_eqn_option) == 'dedt') then
120 44 : s% eps_grav_form_for_energy_eqn = .false.
121 : else
122 0 : call mesa_error(__FILE__,__LINE__,'Invalid choice for energy_eqn_option')
123 : end if
124 :
125 : ! solving structure equations
126 :
127 44 : !$OMP PARALLEL DO PRIVATE(op_err,k) SCHEDULE(dynamic,2)
128 : do k = nzlo, nzhi
129 : op_err = 0
130 : ! hack for composition partials
131 : if (s% fix_d_eos_dxa_partials) then
132 : call fix_d_eos_dxa_partials(s, k, op_err)
133 : if (op_err /= 0) then
134 : if (s% report_ierr) write(*,2) 'ierr from fix_d_eos_dxa_partials', k
135 : ierr = op_err
136 : end if
137 : end if
138 : end do
139 : !$OMP END PARALLEL DO
140 :
141 44 : if (s% use_other_energy_implicit) then
142 0 : call s% other_energy_implicit(s% id, ierr)
143 0 : if (ierr /= 0) then
144 0 : if (len_trim(s% retry_message) == 0) s% retry_message = 'other_energy_implicit failed'
145 0 : if (s% report_ierr) &
146 0 : write(*,*) 'ierr from other_energy_implicit'
147 0 : return
148 : end if
149 : end if
150 :
151 44 : if (s% use_other_momentum_implicit) then
152 0 : call s% other_momentum_implicit(s% id, ierr)
153 0 : if (ierr /= 0) then
154 0 : if (len_trim(s% retry_message) == 0) s% retry_message = 'other_momentum_implicit failed'
155 0 : if (s% report_ierr) &
156 0 : write(*,*) 'ierr from other_momentum_implicit'
157 0 : return
158 : end if
159 : end if
160 :
161 44 : !$OMP PARALLEL DO PRIVATE(op_err,k) SCHEDULE(dynamic,2)
162 : do k = nzlo, nzhi
163 : s% dblk(:,:,k) = 0
164 : s% ublk(:,:,k) = 0
165 : s% lblk(:,:,k) = 0
166 :
167 : op_err = 0
168 : if (do_dlnd_dt) then
169 : call do1_density_eqn(s, k, nvar, op_err)
170 : if (op_err /= 0) then
171 : if (len_trim(s% retry_message) == 0) s% retry_message = 'error in do1_density_eqn'
172 : if (s% report_ierr) write(*,2) 'ierr in do1_density_eqn', k
173 : ierr = op_err
174 : end if
175 : end if
176 : if (k > 1) then ! k=1 is surf P BC
177 : if (do_du_dt) then
178 : call do1_Riemann_momentum_eqn(s, k, nvar, op_err)
179 : if (op_err /= 0) then
180 : if (s% report_ierr) write(*,2) 'ierr in do1_Riemann_momentum_eqn', k
181 : if (len_trim(s% retry_message) == 0) s% retry_message = 'error in do1_Riemann_momentum_eqn'
182 : ierr = op_err
183 : end if
184 : end if
185 : if (do_dv_dt) then
186 : call do1_momentum_eqn(s, k, nvar, op_err)
187 : if (op_err /= 0) then
188 : if (s% report_ierr) write(*,2) 'ierr in do1_momentum_eqn', k
189 : if (len_trim(s% retry_message) == 0) s% retry_message = 'error in do1_momentum_eqn'
190 : ierr = op_err
191 : end if
192 : end if
193 : end if
194 : if (do_dlnR_dt) then
195 : call do1_radius_eqn(s, k, nvar, op_err)
196 : if (op_err /= 0) then
197 : if (s% report_ierr) write(*,2) 'ierr in do1_radius_eqn', k
198 : if (len_trim(s% retry_message) == 0) s% retry_message = 'error in do1_radius_eqn'
199 : ierr = op_err
200 : end if
201 : end if
202 : if (do_alpha_RTI) then
203 : call do1_dalpha_RTI_dt_eqn(s, k, nvar, op_err)
204 : if (op_err /= 0) then
205 : if (s% report_ierr) write(*,2) 'ierr in do1_dalpha_RTI_dt_eqn', k
206 : if (len_trim(s% retry_message) == 0) s% retry_message = 'error in do1_dalpha_RTI_dt_eqn'
207 : ierr = op_err
208 : end if
209 : end if
210 : if (do_w_div_wc) then
211 : call do1_w_div_wc_eqn(s, k, nvar, op_err)
212 : if (op_err /= 0) then
213 : if (s% report_ierr) write(*,2) 'ierr in do1_w_div_wc_eqn', k
214 : if (len_trim(s% retry_message) == 0) s% retry_message = 'error in do1_w_div_wc_eqn'
215 : ierr = op_err
216 : end if
217 : end if
218 : if (do_j_rot) then
219 : call do1_dj_rot_dt_eqn(s, k, nvar, op_err)
220 : if (op_err /= 0) then
221 : if (s% report_ierr) write(*,2) 'ierr in do1_dj_rot_dt_eqn', k
222 : if (len_trim(s% retry_message) == 0) s% retry_message = 'error in do1_dj_rot_dt_eqn'
223 : ierr = op_err
224 : end if
225 : end if
226 :
227 : if (do_dlnE_dt) then
228 : call zero_eps_grav_and_partials(s, k)
229 : call do1_energy_eqn(s, k, do_chem, nvar, op_err)
230 : if (op_err /= 0) then
231 : if (s% report_ierr) write(*,2) 'ierr in do1_energy_eqn', k
232 : if (len_trim(s% retry_message) == 0) s% retry_message = 'error in do1_energy_eqn'
233 : ierr = op_err
234 : end if
235 : end if
236 : if (do_detrb_dt) then
237 : call do1_turbulent_energy_eqn(s, k, nvar, op_err)
238 : if (op_err /= 0) then
239 : if (s% report_ierr) write(*,2) 'ierr in do1_turbulent_energy_eqn', k
240 : if (len_trim(s% retry_message) == 0) s% retry_message = 'error in do1_turbulent_energy_eqn'
241 : ierr = op_err
242 : end if
243 : call do1_rsp2_Hp_eqn(s, k, nvar, op_err)
244 : if (op_err /= 0) then
245 : if (s% report_ierr) write(*,2) 'ierr in do1_rsp2_Hp_eqn', k
246 : if (len_trim(s% retry_message) == 0) s% retry_message = 'error in do1_rsp2_Hp_eqn'
247 : ierr = op_err
248 : end if
249 : end if
250 : if (do_equL) then
251 : if (s% RSP2_flag .and. (k > 1 .or. s% RSP2_use_L_eqn_at_surface)) then
252 : call do1_rsp2_L_eqn(s, k, nvar, op_err)
253 : if (op_err /= 0) then
254 : if (s% report_ierr) write(*,2) 'ierr in do1_rsp2_L_eqn', k
255 : if (len_trim(s% retry_message) == 0) s% retry_message = 'error in do1_rsp2_L_eqn'
256 : ierr = op_err
257 : end if
258 : else if (k > 1) then ! k==1 is done by T_surf BC
259 : call do1_dlnT_dm_eqn(s, k, nvar, op_err)
260 : if (op_err /= 0) then
261 : if (s% report_ierr) write(*,2) 'ierr in do1_dlnT_dm_eqn', k
262 : if (len_trim(s% retry_message) == 0) s% retry_message = 'error in do1_dlnT_dm_eqn'
263 : ierr = op_err
264 : end if
265 : end if
266 : end if
267 :
268 : if (do_chem) then
269 : call do1_chem_eqns(s, k, nvar, op_err)
270 : if (op_err /= 0) then
271 : if (s% report_ierr) write(*,2) 'ierr in do1_chem_eqns', k
272 : if (len_trim(s% retry_message) == 0) s% retry_message = 'error in do1_chem_eqns'
273 : ierr = op_err
274 : end if
275 : end if
276 : end do
277 : !$OMP END PARALLEL DO
278 :
279 44 : if (ierr == 0 .and. nzlo == 1) then
280 44 : call PT_eqns_surf(s, nvar, do_du_dt, do_dv_dt, do_equL, ierr)
281 44 : if (ierr /= 0) then
282 0 : if (s% report_ierr) write(*,2) 'ierr in PT_eqns_surf', ierr
283 0 : if (len_trim(s% retry_message) == 0) s% retry_message = 'error in PT_eqns_surf'
284 : end if
285 : end if
286 :
287 44 : if (ierr /= 0) then
288 0 : if (s% report_ierr) write(*,*) 'ierr in eval_equ_for_solver'
289 0 : return
290 : end if
291 :
292 44 : if (.false. .and. s% model_number == 2) then ! .and. .not. s% doing_relax) then
293 : if (.false.) then
294 : i = s% i_dv_dt
295 : k = 1
296 : do j=1,5
297 : write(*,5) 'dblk(i,j,k) ' // &
298 : trim(s% nameofequ(i)) // ' ' // trim(s% nameofvar(j)), &
299 : s% solver_iter, i, j, k, s% dblk(i,j,k)
300 : end do
301 : end if
302 : !write(*,*) 'call show_matrix'
303 : !call show_matrix(s, s% dblk(1:s% nvar_hydro,1:s% nvar_hydro,1), s% nvar_hydro)
304 : call dump_equ ! debugging
305 : call mesa_error(__FILE__,__LINE__,'after dump_equ')
306 : end if
307 :
308 :
309 : contains
310 :
311 : subroutine dump_equ
312 : integer :: k, j
313 : include 'formats'
314 : do k=1,s% nz
315 : do j=1,nvar
316 : write(*,3) 'equ ' // trim(s% nameofequ(j)), &
317 : k, s% solver_iter, s% equ(j, k)
318 : end do
319 : write(*,3) 'eps_mdot', k, s% solver_iter, s% eps_mdot(k)
320 : write(*,3) 'energy_start', k, s% solver_iter, s% energy_start(k)
321 : write(*,3) 'energy', k, s% solver_iter, s% energy(k)
322 : write(*,3) 'L', k, s% solver_iter, s% L(k)
323 : write(*,3) 'L', k+1, s% solver_iter, s% L(k+1)
324 : write(*,'(A)')
325 : !if (k == 6) exit
326 : end do
327 44 : end subroutine dump_equ
328 :
329 44 : subroutine unpack
330 :
331 : include 'formats'
332 :
333 44 : nz = s% nz
334 44 : species = s% species
335 44 : nvar_hydro = s% nvar_hydro
336 :
337 44 : lnT => s% lnT
338 44 : lnR => s% lnR
339 44 : L => s% L
340 44 : lnP => s% lnPeos
341 44 : energy => s% energy
342 :
343 44 : i_dv_dt = s% i_dv_dt
344 44 : i_du_dt = s% i_du_dt
345 44 : i_equL = s% i_equL
346 44 : i_dlnd_dt = s% i_dlnd_dt
347 44 : i_dlnE_dt = s% i_dlnE_dt
348 44 : i_dlnR_dt = s% i_dlnR_dt
349 44 : i_dalpha_RTI_dt = s% i_dalpha_RTI_dt
350 44 : i_equ_w_div_wc = s% i_equ_w_div_wc
351 44 : i_dj_rot_dt = s% i_dj_rot_dt
352 44 : i_detrb_dt = s% i_detrb_dt
353 :
354 44 : i_lnd = s% i_lnd
355 44 : i_lnT = s% i_lnT
356 44 : i_lnR = s% i_lnR
357 44 : i_lum = s% i_lum
358 44 : i_v = s% i_v
359 44 : i_u = s% i_u
360 44 : i_alpha_RTI = s% i_alpha_RTI
361 44 : i_w_div_wc = s% i_w_div_wc
362 44 : i_j_rot = s% i_j_rot
363 :
364 44 : i_xh1 = s% nvar_hydro+s% net_iso(ih1)
365 44 : i_xhe4 = s% nvar_hydro+s% net_iso(ihe4)
366 :
367 44 : L_phot_old = s% L_phot_old
368 44 : v_flag = s% v_flag
369 44 : u_flag = s% u_flag
370 :
371 44 : end subroutine unpack
372 :
373 52348 : subroutine fix_d_eos_dxa_partials(s, k, ierr)
374 :
375 : ! revise composition partials
376 : ! subroutine can be removed when EOS fully provides composition partials
377 :
378 : use eos_def, only: num_eos_basic_results, num_eos_d_dxa_results, i_lnE, i_lnPgas
379 : use eos_support, only: get_eos
380 : use star_utils, only: lookup_nameofvar
381 : type (star_info), pointer :: s
382 : integer, intent(in) :: k
383 : integer, intent(out) :: ierr
384 :
385 : integer :: j
386 :
387 : logical, parameter :: debug = .false.
388 :
389 : ! these vars are for faking composition derivatives
390 : real(dp), dimension(num_eos_basic_results) :: &
391 4187840 : res, dres_dlnd, dres_dlnT
392 1256352 : real(dp) :: dres_dxa(num_eos_d_dxa_results, s% species)
393 471132 : real(dp) :: dxa
394 471132 : real(dp) :: xa_start_1(s% species)
395 471132 : real(dp) :: frac_without_dxa
396 471132 : real(dp) :: lnE_with_xa_start, lnPgas_with_xa_start
397 :
398 : integer :: i_var, i_var_sink
399 :
400 : real(dp), parameter :: dxa_threshold = 1d-4
401 :
402 : logical, parameter :: checking = .true.
403 :
404 : include 'formats'
405 :
406 52348 : ierr = 0
407 :
408 : ! some EOSes have composition partials and some do not
409 : ! those currently without dx partials are PC & Skye & ideal
410 52348 : frac_without_dxa = s% eos_frac_PC(k) + s% eos_frac_Skye(k) + s% eos_frac_ideal(k)
411 :
412 : if (debug .and. k == s% solver_test_partials_k) then
413 : write(*,2) 's% eos_frac_PC(k)', k, s% eos_frac_PC(k)
414 : write(*,2) 's% eos_frac_Skye(k)', k, s% eos_frac_Skye(k)
415 : write(*,2) 's% eos_frac_ideal(k)', k, s% eos_frac_ideal(k)
416 : write(*,2) 'frac_without_dxa', k, frac_without_dxa
417 : end if
418 :
419 52348 : if (k == s% solver_test_partials_k .and. s% solver_iter == s% solver_test_partials_iter_number) then
420 0 : i_var = lookup_nameofvar(s, s% solver_test_partials_var_name)
421 0 : if (i_var > s% nvar_hydro) then
422 0 : i_var_sink = lookup_nameofvar(s, s% solver_test_partials_sink_name)
423 : end if
424 : end if
425 :
426 : ! if we're on an EOS where there aren't composition partials,
427 : ! approximate derivatives with finite differences
428 52348 : if (frac_without_dxa > 0) then
429 :
430 43704 : do j=1, s% species
431 38848 : dxa = s% xa(j,k) - s% xa_start(j,k)
432 :
433 : if (debug .and. k == s% solver_test_partials_k .and. &
434 : s% solver_iter == s% solver_test_partials_iter_number) &
435 : write(*,2) 'dxa', j, dxa
436 :
437 43704 : if (abs(dxa) >= dxa_threshold) then
438 :
439 : ! first, get eos with xa_start
440 :
441 : call get_eos( &
442 : s, k, s% xa_start(:,k), &
443 : s% rho(k), s% lnd(k)/ln10, s% T(k), s% lnT(k)/ln10, &
444 0 : res, dres_dlnd, dres_dlnT, dres_dxa, ierr)
445 0 : if (ierr /= 0) then
446 0 : if (s% report_ierr) write(*,2) 'failed in get_eos with xa_start', k
447 0 : return
448 : end if
449 :
450 0 : lnE_with_xa_start = res(i_lnE)
451 0 : lnPgas_with_xa_start = res(i_lnPgas)
452 :
453 : ! now, get eos with 1 iso perturbed
454 :
455 0 : xa_start_1 = s% xa_start(:,k)
456 0 : xa_start_1(j) = s% xa_start(j,k) + dxa
457 :
458 : call get_eos( &
459 : s, k, xa_start_1, &
460 : s% Rho(k), s% lnd(k)/ln10, s% T(k), s% lnT(k)/ln10, &
461 0 : res, dres_dlnd, dres_dlnT, dres_dxa, ierr)
462 0 : if (is_bad(res(i_lnPgas))) ierr = -1
463 0 : if (ierr /= 0) then
464 :
465 : ! punt silently for now
466 0 : s% dlnE_dxa_for_partials(:,k) = 0d0
467 0 : s% dlnPeos_dxa_for_partials(:,k) = 0d0
468 0 : ierr = 0
469 0 : return
470 :
471 : if (s% report_ierr) write(*,2) 'failed in get_eos with xa_start_1', k
472 : return
473 : end if
474 :
475 : ! fix up derivatives
476 :
477 : if (debug .and. k == s% solver_test_partials_k) & ! .and. s% solver_iter == s% solver_test_partials_iter_number) &
478 : write(*,2) 'res(i_lnE) - lnE_with_xa_start', j, res(i_lnE) - lnE_with_xa_start
479 :
480 : s% dlnE_dxa_for_partials(j,k) = dres_dxa(i_lnE, j) + &
481 0 : frac_without_dxa * (res(i_lnE) - lnE_with_xa_start) / dxa
482 0 : if (checking) call check_dxa(j,k,s% dlnE_dxa_for_partials(j,k),'dlnE_dxa_for_partials')
483 :
484 : s% dlnPeos_dxa_for_partials(j,k) = s% Pgas(k)*dres_dxa(i_lnPgas,j)/s% Peos(k) + &
485 0 : frac_without_dxa * (s% Pgas(k)/s% Peos(k)) * (res(i_lnPgas) - lnPgas_with_xa_start) / dxa
486 : if (.false. .and. s% model_number == 1100 .and. k == 151 .and. &
487 : j == 1 .and. is_bad(s% dlnPeos_dxa_for_partials(j,k))) then
488 : write(*,2) 's% Pgas(k)', k, s% Pgas(k)
489 : write(*,2) 'dres_dxa(i_lnPgas,j)', k, dres_dxa(i_lnPgas,j)
490 : write(*,2) 's% Peos(k)', k, s% Peos(k)
491 : write(*,2) 'frac_without_dxa', k, frac_without_dxa
492 : write(*,2) 'res(i_lnPgas)', k, res(i_lnPgas)
493 : write(*,2) 'lnPgas_with_xa_start', k, lnPgas_with_xa_start
494 : write(*,2) 'dxa', k, dxa
495 : write(*,2) 'dxa_threshold', k, dxa_threshold
496 : write(*,2) 's% xa_start(j,k)', j, s% xa_start(j,k)
497 : write(*,2) 'xa_start_1(j)', j, xa_start_1(j)
498 : write(*,2) 's% eos_frac_PC(k)', k, s% eos_frac_PC(k)
499 : write(*,2) 's% eos_frac_Skye(k)', k, s% eos_frac_Skye(k)
500 : end if
501 0 : if (checking) call check_dxa(j,k,s% dlnPeos_dxa_for_partials(j,k),'dlnPeos_dxa_for_partials')
502 :
503 : else
504 :
505 38848 : if (k == s% solver_test_partials_k .and. s% solver_iter == s% solver_test_partials_iter_number) then
506 0 : if (i_var_sink > 0 .and. i_var > s% nvar_hydro) then
507 0 : if (dxa < dxa_threshold) then
508 0 : if (j == i_var - s% nvar_hydro) then
509 0 : write(*,*) 'fix_d_eos_dxa_partials: skipping dxa derivative fix for ', &
510 0 : trim (s% solver_test_partials_var_name), &
511 0 : ' (dxa < dxa_threshold): ', abs(dxa), ' < ', dxa_threshold
512 : end if
513 0 : if (j == i_var_sink - s% nvar_hydro) then
514 0 : write(*,*) 'fix_d_eos_dxa_partials: skipping dxa derivative fix for ', &
515 0 : trim (s% solver_test_partials_sink_name), &
516 0 : ' (dxa < dxa_threshold): ', abs(dxa), ' < ', dxa_threshold
517 : end if
518 : end if
519 : end if
520 : end if
521 :
522 : end if
523 : end do
524 :
525 : end if
526 :
527 52348 : end subroutine fix_d_eos_dxa_partials
528 :
529 0 : subroutine check_dxa(j, k, dxa, str)
530 : integer, intent(in) :: j, k
531 : real(dp), intent(in) :: dxa
532 : character (len=*), intent(in) :: str
533 : include 'formats'
534 0 : if (is_bad(dxa)) then
535 0 : !$omp critical (eval_equ_for_solver_crit1)
536 0 : ierr = -1
537 0 : if (s% report_ierr) then
538 0 : write(*,3) 'eval_equ_for_solver: bad ' // trim(str), j, k, dxa
539 : end if
540 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'eval_equ_for_solver')
541 : !$omp end critical (eval_equ_for_solver_crit1)
542 0 : return
543 : end if
544 52348 : end subroutine check_dxa
545 :
546 :
547 : end subroutine eval_equ_for_solver
548 :
549 :
550 52348 : subroutine do1_density_eqn(s, k, nvar, ierr)
551 : use star_utils, only: save_eqn_residual_info
552 : type (star_info), pointer :: s
553 : integer, intent(in) :: k, nvar
554 : integer, intent(out) :: ierr
555 :
556 : type(auto_diff_real_star_order1) :: &
557 : rho, dr3, r_p1, rp13, lnR_expected, lnR_actual, resid_ad
558 : logical :: test_partials
559 : include 'formats'
560 :
561 : !test_partials = (k == s% solver_test_partials_k)
562 52348 : test_partials = .false.
563 :
564 52348 : ierr = 0
565 52348 : rho = wrap_d_00(s,k)
566 52348 : dr3 = (s% dm(k)/rho)/(4d0*pi/3d0)
567 52348 : r_p1 = wrap_r_p1(s,k)
568 52348 : rp13 = pow3(r_p1)
569 :
570 52348 : lnR_expected = log(rp13 + dr3)/3d0
571 52348 : lnR_actual = wrap_lnR_00(s,k)
572 :
573 52348 : resid_ad = lnR_actual - lnR_expected
574 52348 : s% equ(s% i_dlnd_dt, k) = resid_ad%val
575 :
576 : if (test_partials) then
577 : s% solver_test_partials_val = 0
578 : end if
579 :
580 : call save_eqn_residual_info( &
581 52348 : s, k, nvar, s% i_dlnd_dt, resid_ad, 'do1_density_eqn', ierr)
582 :
583 : if (test_partials) then
584 : s% solver_test_partials_var = 0
585 : s% solver_test_partials_dval_dx = 0
586 : write(*,*) 'do1_density_eqn', s% solver_test_partials_var
587 : end if
588 :
589 52348 : end subroutine do1_density_eqn
590 :
591 :
592 0 : subroutine do1_w_div_wc_eqn(s, k, nvar, ierr)
593 52348 : use hydro_rotation
594 : use star_utils, only: save_eqn_residual_info
595 : type (star_info), pointer :: s
596 : integer, intent(in) :: k, nvar
597 : integer, intent(out) :: ierr
598 : integer :: i_equ_w_div_wc, i_w_div_wc
599 : real(dp) :: wwc, wwc4, wwc6, wwc_log_term, dimless_rphi, dimless_rphi_given_wwc, w1, w2, jr_lim1, jr_lim2
600 : type(auto_diff_real_star_order1) :: &
601 : w_d_wc00, w4_d_wc00, w6_d_wc00, r00, w_log_term_d_wc00, jrot00, resid_ad, A_ad, C_ad, &
602 : jrot_ratio, sigmoid_jrot_ratio
603 : logical :: test_partials
604 : include 'formats'
605 :
606 0 : ierr = 0
607 :
608 : !test_partials = (k == s% solver_test_partials_k-1)
609 0 : test_partials = .false.
610 :
611 0 : i_equ_w_div_wc = s% i_equ_w_div_wc
612 0 : i_w_div_wc = s% i_w_div_wc
613 :
614 0 : wwc = s% w_div_wcrit_max
615 0 : wwc4 = pow4(wwc)
616 0 : wwc6 = pow6(wwc)
617 0 : wwc_log_term = log(1d0 - pow(wwc, log_term_power))
618 0 : jr_lim1 = two_thirds * wwc * C(pow2(wwc), wwc4, wwc6, wwc_log_term) / A(wwc4, wwc6, wwc_log_term)
619 :
620 0 : wwc = s% w_div_wcrit_max2
621 0 : wwc4 = pow4(wwc)
622 0 : wwc6 = pow6(wwc)
623 0 : wwc_log_term = log(1d0 - pow(wwc, log_term_power))
624 0 : jr_lim2 = two_thirds * wwc * C(pow2(wwc), wwc4, wwc6, wwc_log_term) / A(wwc4, wwc6, wwc_log_term)
625 :
626 0 : w_d_wc00 = wrap_w_div_wc_00(s, k)
627 0 : w4_d_wc00 = pow4(w_d_wc00)
628 0 : w6_d_wc00 = pow6(w_d_wc00)
629 0 : w_log_term_d_wc00 = log(1d0 - pow(w_d_wc00, log_term_power))
630 0 : A_ad = 1d0 + 0.3293d0 * w4_d_wc00 - 0.4926d0 * w6_d_wc00 - 0.5560d0 * w_log_term_d_wc00
631 : C_ad = 1d0 + 17d0 / 60d0 * pow2(w_d_wc00) + 0.4010d0 * w4_d_wc00 - 0.8606d0 * w6_d_wc00 &
632 0 : - 0.9305d0 * w_log_term_d_wc00
633 :
634 0 : r00 = wrap_r_00(s, k)
635 0 : if (s% j_rot_flag) then
636 0 : jrot00 = wrap_jrot_00(s, k)
637 0 : jrot_ratio = jrot00 / sqrt(s% cgrav(k) * s% m(k) * r00)
638 : else
639 0 : jrot_ratio = s% j_rot(k) / sqrt(s% cgrav(k) * s% m(k) * r00)
640 : end if
641 0 : if (abs(jrot_ratio% val) > jr_lim1) then
642 : sigmoid_jrot_ratio = 2d0 * (jr_lim2-jr_lim1) &
643 : / (1d0 + exp(-2d0 * (abs(jrot_ratio) - jr_lim1) / (jr_lim2 - jr_lim1))) &
644 0 : - jr_lim2 + 2 * jr_lim1
645 0 : if (jrot_ratio% val < 0d0) then
646 0 : sigmoid_jrot_ratio = -sigmoid_jrot_ratio
647 : end if
648 0 : resid_ad = (sigmoid_jrot_ratio - two_thirds * w_d_wc00 * C_ad / A_ad)
649 : else
650 0 : resid_ad = (jrot_ratio - two_thirds * w_d_wc00 * C_ad / A_ad)
651 : end if
652 :
653 0 : s% equ(i_equ_w_div_wc, k) = resid_ad% val
654 :
655 : call save_eqn_residual_info( &
656 0 : s, k, nvar, i_equ_w_div_wc, resid_ad, 'do1_w_div_wc_eqn', ierr)
657 :
658 : if (test_partials) then
659 : s% solver_test_partials_var = s% i_w_div_wc
660 : s% solver_test_partials_dval_dx = 1d0
661 : write(*,*) 'do1_w_div_wc_eqn', s% solver_test_partials_var
662 : end if
663 :
664 0 : end subroutine do1_w_div_wc_eqn
665 :
666 :
667 0 : subroutine do1_dj_rot_dt_eqn(s, k, nvar, ierr)
668 0 : use hydro_rotation
669 : use star_utils, only: save_eqn_residual_info
670 : type (star_info), pointer :: s
671 : integer, intent(in) :: k, nvar
672 : integer, intent(out) :: ierr
673 : integer :: i_dj_rot_dt, i_j_rot
674 : real(dp) :: scale
675 : type(auto_diff_real_star_order1) :: resid_ad, jrot00, djrot_dt, F00, Fm1, flux_diff
676 : logical :: test_partials
677 : include 'formats'
678 :
679 0 : ierr = 0
680 :
681 : !test_partials = (k == s% solver_test_partials_k)
682 0 : test_partials = .false.
683 :
684 : ! Need to find a good way to scale the equation
685 : scale = 1d6*max(1d2*sqrt(s% cgrav(k)*s% m(k)*s%r_start(k))/s%dt,&
686 0 : s% total_abs_angular_momentum/(s% dm(k)*s% dt*s% nz))
687 :
688 0 : i_dj_rot_dt = s% i_dj_rot_dt
689 0 : i_j_rot = s% i_j_rot
690 :
691 0 : jrot00 = wrap_jrot_00(s, k)
692 0 : F00 = s% j_flux(k)
693 0 : if (k > 1) then
694 0 : Fm1 = shift_m1(s% j_flux(k-1))
695 : else
696 0 : Fm1 = 0d0
697 : end if
698 :
699 0 : djrot_dt = (jrot00-s% j_rot_start(k))/s% dt
700 0 : flux_diff = (Fm1-F00)/s% dm_bar(k)
701 :
702 0 : resid_ad = (djrot_dt+flux_diff-s% extra_jdot(k))/scale
703 :
704 0 : s% equ(i_dj_rot_dt, k) = resid_ad% val
705 :
706 : call save_eqn_residual_info( &
707 0 : s, k, nvar, i_dj_rot_dt, resid_ad, 'do1_dj_rot_dt_eqn', ierr)
708 :
709 : !if (test_partials) then
710 : ! s% solver_test_partials_val = s% i_rot(k)
711 : !end if
712 :
713 : !if (test_partials) then
714 : ! s% solver_test_partials_var = s% i_lnR
715 : ! s% solver_test_partials_dval_dx = s% di_rot_dlnR(k)
716 : ! write(*,*) 'do1_w_div_wc_eqn', s% solver_test_partials_var
717 : !end if
718 :
719 0 : end subroutine do1_dj_rot_dt_eqn
720 :
721 :
722 44 : subroutine PT_eqns_surf(s, nvar, do_du_dt, do_dv_dt, do_equL, ierr)
723 :
724 0 : use star_utils, only: save_eqn_residual_info
725 : use eos_lib, only: Radiation_Pressure
726 :
727 : type (star_info), pointer :: s
728 : integer, intent(in) :: nvar
729 : logical, intent(in) :: do_du_dt, do_dv_dt, do_equL
730 : integer, intent(out) :: ierr
731 :
732 : type(auto_diff_real_star_order1) :: &
733 : P_bc_ad, T_bc_ad, lnP_bc_ad, lnT_bc_ad, resid_ad
734 : integer :: i_P_eqn
735 : logical :: offset_T_to_cell_center, offset_P_to_cell_center, &
736 : need_P_surf, need_T_surf, test_partials
737 :
738 : include 'formats'
739 :
740 : !test_partials = (s% solver_iter == s% solver_test_partials_iter_number)
741 44 : test_partials = .false.
742 44 : ierr = 0
743 44 : if (s% u_flag) then
744 0 : i_P_eqn = s% i_du_dt
745 : else ! use this even if not v_flag
746 44 : i_P_eqn = s% i_dv_dt
747 : end if
748 :
749 :
750 44 : if(s% drag_coefficient > 0) then
751 : ! We dont call expected_non_HSE_term with k==1 unless we call set_momentum_BC
752 : ! so lets initilize this to zero, then if we dont call set_momentum_BC we have a
753 : ! sensible value here.
754 0 : s% dvdt_drag(1) = 0
755 : end if
756 :
757 44 : need_P_surf = .false.
758 44 : if (s% use_compression_outer_BC) then
759 0 : call set_compression_BC(ierr)
760 44 : else if (s% use_zero_Pgas_outer_BC) then
761 0 : P_bc_ad = Radiation_Pressure(s% T_start(1))
762 0 : call set_momentum_BC(ierr)
763 44 : else if (s% use_fixed_Psurf_outer_BC) then
764 : P_bc_ad = s% fixed_Psurf
765 0 : call set_momentum_BC(ierr)
766 44 : else if (s% use_fixed_vsurf_outer_BC) then
767 0 : call set_fixed_vsurf_outer_BC(ierr)
768 : else
769 44 : need_P_surf = .true.
770 : end if
771 44 : if (ierr /= 0) return
772 :
773 44 : need_T_surf = .false.
774 44 : if ((.not. do_equL) .or. &
775 : (s% RSP2_flag .and. s% RSP2_use_L_eqn_at_surface)) then
776 : ! no Tsurf BC
777 : else
778 44 : need_T_surf = .true.
779 : end if
780 : if (ierr /= 0) return
781 :
782 44 : offset_P_to_cell_center = .not. s% use_momentum_outer_BC
783 :
784 44 : offset_T_to_cell_center = .true.
785 44 : if (s% use_other_surface_PT .or. s% RSP2_flag) &
786 0 : offset_T_to_cell_center = .false.
787 :
788 44 : call get_PT_bc_ad(ierr)
789 44 : if (ierr /= 0) return
790 :
791 44 : if (need_P_surf) then
792 44 : if (s% use_momentum_outer_BC) then
793 0 : call set_momentum_BC(ierr)
794 : else
795 44 : call set_Psurf_BC(ierr)
796 : end if
797 44 : if (ierr /= 0) return
798 : end if
799 :
800 44 : if (need_T_surf) then
801 44 : call set_Tsurf_BC(ierr)
802 44 : if (ierr /= 0) return
803 : end if
804 :
805 : if (test_partials) then
806 : s% solver_test_partials_val = 0
807 : end if
808 :
809 44 : if (test_partials) then
810 : s% solver_test_partials_var = 0
811 : s% solver_test_partials_dval_dx = 0
812 : write(*,*) 'PT_eqns_surf', s% solver_test_partials_var
813 : end if
814 :
815 : contains
816 :
817 44 : subroutine get_PT_bc_ad(ierr) ! set P_bc_ad and T_bc_ad
818 44 : use hydro_vars, only: set_Teff_info_for_eqns
819 : use chem_def
820 : use atm_def
821 : integer, intent(out) :: ierr
822 : real(dp) :: r, L, Teff, &
823 : lnT_surf, dlnTsurf_dL, dlnTsurf_dlnR, dlnTsurf_dlnM, dlnTsurf_dlnkap, &
824 : lnP_surf, dlnPsurf_dL, dlnPsurf_dlnR, dlnPsurf_dlnM, dlnPsurf_dlnkap
825 : real(dp) :: &
826 44 : dlnT_bc_dlnd, dlnT_bc_dlnT, dlnT_bc_dlnR, &
827 44 : dlnT_bc_dL, dlnP_bc_dlnd, dlnP_bc_dlnT, dlnP_bc_dL, dlnP_bc_dlnR, &
828 44 : dlnkap_dlnd, dlnkap_dlnT, dPinv_dlnd, dPinv_dlnT, dP0, dT0, &
829 44 : P_surf, T_surf, dlnP_bc_dlnPsurf, &
830 44 : dlnT_bc_dlnTsurf, P_bc, T_bc, lnT_bc, lnP_bc, &
831 44 : dP0_dlnR, dT0_dlnR, dT0_dlnT, dT0_dlnd, dT0_dL, dlnP_bc_dP0, dlnT_bc_dT0, &
832 44 : d_gradT_dlnR, d_gradT_dlnT00, d_gradT_dlnd00, d_gradT_dL, &
833 44 : dlnR00, dlnT00, dlnd00
834 : logical, parameter :: skip_partials = .false.
835 : include 'formats'
836 : ierr = 0
837 :
838 : call set_Teff_info_for_eqns(s, skip_partials, &
839 : need_P_surf, need_T_surf, r, L, Teff, &
840 : lnT_surf, dlnTsurf_dL, dlnTsurf_dlnR, dlnTsurf_dlnM, dlnTsurf_dlnkap, &
841 : lnP_surf, dlnPsurf_dL, dlnPsurf_dlnR, dlnPsurf_dlnM, dlnPsurf_dlnkap, &
842 44 : ierr)
843 44 : if (ierr /= 0) then
844 0 : if (s% report_ierr) then
845 0 : write(*,*) 'PT_eqns_surf: ierr from set_Teff_info_for_eqns'
846 : end if
847 : return
848 : end if
849 :
850 : ! P_surf and T_surf are at outer boundary of cell 1
851 44 : P_surf = exp(lnP_surf)
852 44 : T_surf = exp(lnT_surf)
853 :
854 44 : s% P_surf = P_surf
855 44 : s% T_surf = T_surf
856 :
857 44 : dP0 = 0
858 44 : dT0 = 0
859 44 : if (offset_P_to_cell_center) &
860 44 : dP0 = s% cgrav(1)*s% m_grav(1)*s% dm(1)/(8*pi*pow4(r))
861 44 : if (offset_T_to_cell_center) &
862 44 : dT0 = dP0*s% gradT(1)*s% T(1)/s% Peos(1)
863 :
864 44 : P_bc = P_surf + dP0
865 44 : T_bc = T_surf + dT0
866 :
867 44 : lnP_bc = log(P_bc)
868 44 : lnT_bc = log(T_bc)
869 :
870 44 : if (is_bad(P_bc)) then
871 0 : write(*,1) 'lnP_bc', lnP_bc
872 0 : write(*,1) 'P_bc', P_bc
873 0 : write(*,1) 'P_surf', P_surf
874 0 : write(*,1) 'dP0', dP0
875 0 : write(*,1) 'lnP_surf', lnP_surf
876 0 : write(*,1) 'r', r
877 0 : call mesa_error(__FILE__,__LINE__,'P bc')
878 : end if
879 :
880 44 : if (is_bad(T_bc)) then
881 0 : write(*,1) 'lnT_bc', lnT_bc
882 0 : write(*,1) 'T_bc', T_bc
883 0 : write(*,1) 'T_surf', T_surf
884 0 : write(*,1) 'dP0', dP0
885 0 : write(*,1) 'lnT_surf', lnT_surf
886 0 : call mesa_error(__FILE__,__LINE__,'T bc')
887 : end if
888 :
889 44 : dP0_dlnR = 0
890 44 : if (offset_P_to_cell_center) then ! include partials of dP0
891 44 : dP0_dlnR = -4*dP0
892 : end if
893 :
894 44 : dT0_dlnR = 0
895 44 : dT0_dlnT = 0
896 44 : dT0_dlnd = 0
897 44 : dT0_dL = 0
898 44 : if (offset_T_to_cell_center) then ! include partials of dT0
899 44 : d_gradT_dlnR = s% gradT_ad(1)%d1Array(i_lnR_00)
900 44 : d_gradT_dlnT00 = s% gradT_ad(1)%d1Array(i_lnT_00)
901 44 : d_gradT_dlnd00 = s% gradT_ad(1)%d1Array(i_lnd_00)
902 44 : d_gradT_dL = s% gradT_ad(1)%d1Array(i_L_00)
903 44 : dT0_dlnR = -4*dT0 + dP0*d_gradT_dlnR*s% T(1)/s% Peos(1)
904 44 : dPinv_dlnT = -s% chiT_for_partials(1)/s% Peos(1)
905 : dT0_dlnT = &
906 : dT0 + &
907 : dP0*d_gradT_dlnT00*s% T(1)/s% Peos(1) + &
908 44 : dP0*s% gradT(1)*s% T(1)*dPinv_dlnT
909 44 : dPinv_dlnd = -s% chiRho_for_partials(1)/s% Peos(1)
910 : dT0_dlnd = &
911 : dP0*d_gradT_dlnd00*s% T(1)/s% Peos(1) + &
912 44 : dP0*s% gradT(1)*s% T(1)*dPinv_dlnd
913 44 : dT0_dL = dP0*d_gradT_dL*s% T(1)/s% Peos(1)
914 : end if
915 :
916 44 : dlnP_bc_dP0 = 1/P_bc
917 44 : dlnT_bc_dT0 = 1/T_bc
918 :
919 44 : dlnP_bc_dlnPsurf = P_surf/P_bc
920 44 : dlnT_bc_dlnTsurf = T_surf/T_bc
921 :
922 44 : dlnkap_dlnd = s% d_opacity_dlnd(1)/s% opacity(1)
923 44 : dlnkap_dlnT = s% d_opacity_dlnT(1)/s% opacity(1)
924 :
925 44 : dlnP_bc_dlnd = dlnP_bc_dlnPsurf*dlnPsurf_dlnkap*dlnkap_dlnd
926 44 : dlnP_bc_dlnT = dlnP_bc_dlnPsurf*dlnPsurf_dlnkap*dlnkap_dlnT
927 44 : dlnP_bc_dL = dlnP_bc_dlnPsurf*dlnPsurf_dL
928 44 : dlnP_bc_dlnR = dlnP_bc_dlnPsurf*dlnPsurf_dlnR + dlnP_bc_dP0*dP0_dlnR
929 :
930 44 : dlnR00 = P_bc*dlnP_bc_dlnR
931 44 : dlnT00 = P_bc*dlnP_bc_dlnT
932 44 : dlnd00 = P_bc*dlnP_bc_dlnd
933 : call wrap(P_bc_ad, P_bc, &
934 : 0d0, dlnd00, 0d0, &
935 : 0d0, dlnT00, 0d0, &
936 : 0d0, 0d0, 0d0, &
937 : 0d0, dlnR00, 0d0, &
938 : 0d0, 0d0, 0d0, &
939 : 0d0, P_bc*dlnP_bc_dL, 0d0, &
940 : 0d0, 0d0, 0d0, &
941 : 0d0, 0d0, 0d0, &
942 : 0d0, 0d0, 0d0, &
943 : 0d0, 0d0, 0d0, &
944 44 : 0d0, 0d0, 0d0)
945 44 : dlnR00 = dlnP_bc_dlnR
946 44 : dlnT00 = dlnP_bc_dlnT
947 44 : dlnd00 = dlnP_bc_dlnd
948 : call wrap(lnP_bc_ad, lnP_bc, &
949 : 0d0, dlnd00, 0d0, &
950 : 0d0, dlnT00, 0d0, &
951 : 0d0, 0d0, 0d0, &
952 : 0d0, dlnR00, 0d0, &
953 : 0d0, 0d0, 0d0, &
954 : 0d0, dlnP_bc_dL, 0d0, &
955 : 0d0, 0d0, 0d0, &
956 : 0d0, 0d0, 0d0, &
957 : 0d0, 0d0, 0d0, &
958 : 0d0, 0d0, 0d0, &
959 44 : 0d0, 0d0, 0d0)
960 :
961 : dlnT_bc_dlnT = dlnT_bc_dlnTsurf*dlnTsurf_dlnkap*dlnkap_dlnT &
962 44 : + dlnT_bc_dT0*dT0_dlnT
963 : dlnT_bc_dlnd = dlnT_bc_dlnTsurf*dlnTsurf_dlnkap*dlnkap_dlnd &
964 44 : + dlnT_bc_dT0*dT0_dlnd
965 44 : dlnT_bc_dL = dlnT_bc_dlnTsurf*dlnTsurf_dL + dlnT_bc_dT0*dT0_dL
966 44 : dlnT_bc_dlnR = dlnT_bc_dlnTsurf*dlnTsurf_dlnR + dlnT_bc_dT0*dT0_dlnR
967 :
968 44 : dlnR00 = T_bc*dlnT_bc_dlnR
969 44 : dlnT00 = T_bc*dlnT_bc_dlnT
970 44 : dlnd00 = T_bc*dlnT_bc_dlnd
971 : call wrap(T_bc_ad, T_bc, &
972 : 0d0, dlnd00, 0d0, &
973 : 0d0, dlnT00, 0d0, &
974 : 0d0, 0d0, 0d0, &
975 : 0d0, dlnR00, 0d0, &
976 : 0d0, 0d0, 0d0, &
977 : 0d0, T_bc*dlnT_bc_dL, 0d0, &
978 : 0d0, 0d0, 0d0, &
979 : 0d0, 0d0, 0d0, &
980 : 0d0, 0d0, 0d0, &
981 : 0d0, 0d0, 0d0, &
982 44 : 0d0, 0d0, 0d0)
983 44 : dlnR00 = dlnT_bc_dlnR
984 44 : dlnT00 = dlnT_bc_dlnT
985 44 : dlnd00 = dlnT_bc_dlnd
986 : call wrap(lnT_bc_ad, lnT_bc, &
987 : 0d0, dlnd00, 0d0, &
988 : 0d0, dlnT00, 0d0, &
989 : 0d0, 0d0, 0d0, &
990 : 0d0, dlnR00, 0d0, &
991 : 0d0, 0d0, 0d0, &
992 : 0d0, dlnT_bc_dL, 0d0, &
993 : 0d0, 0d0, 0d0, &
994 : 0d0, 0d0, 0d0, &
995 : 0d0, 0d0, 0d0, &
996 : 0d0, 0d0, 0d0, &
997 44 : 0d0, 0d0, 0d0)
998 :
999 44 : end subroutine get_PT_bc_ad
1000 :
1001 :
1002 44 : subroutine set_Tsurf_BC(ierr)
1003 : integer, intent(out) :: ierr
1004 : logical :: test_partials
1005 : type(auto_diff_real_star_order1) :: &
1006 : lnT1_ad, dT4_dm, T4_p1, T4_surf, T4_00_actual, T4_00_expected
1007 44 : real(dp) :: residual
1008 : include 'formats'
1009 : !test_partials = (1 == s% solver_test_partials_k)
1010 44 : test_partials = .false.
1011 44 : ierr = 0
1012 44 : if (s% RSP2_flag) then ! interpolate lnT by mass
1013 0 : T4_p1 = pow4(wrap_T_p1(s,1))
1014 0 : T4_surf = pow4(T_bc_ad)
1015 0 : dT4_dm = (T4_surf - T4_p1)/(s% dm(1) + 0.5d0*s% dm(2))
1016 0 : T4_00_expected = T4_surf - 0.5d0*s% dm(1)*dT4_dm
1017 0 : T4_00_actual = pow4(wrap_T_00(s,1))
1018 0 : resid_ad = T4_00_expected/T4_00_actual - 1d0
1019 : else
1020 44 : lnT1_ad = wrap_lnT_00(s,1)
1021 44 : resid_ad = lnT_bc_ad/lnT1_ad - 1d0
1022 : end if
1023 44 : residual = resid_ad%val
1024 44 : s% equ(s% i_equL, 1) = residual
1025 44 : if (is_bad(residual)) then
1026 0 : write(*,1) 'set_Tsurf_BC residual', residual
1027 0 : write(*,1) 'lnT1_ad%val', lnT1_ad%val
1028 0 : write(*,1) 'lnT_bc_ad%val', lnT_bc_ad%val
1029 0 : call mesa_error(__FILE__,__LINE__,'set_Tsurf_BC')
1030 : end if
1031 : if (test_partials) then
1032 : s% solver_test_partials_val = 0
1033 : end if
1034 : call save_eqn_residual_info( &
1035 44 : s, 1, nvar, s% i_equL, resid_ad, 'set_Tsurf_BC', ierr)
1036 : if (test_partials) then
1037 : s% solver_test_partials_var = 0
1038 : s% solver_test_partials_dval_dx = 0
1039 : write(*,*) 'set_Tsurf_BC', s% solver_test_partials_var
1040 : end if
1041 44 : end subroutine set_Tsurf_BC
1042 :
1043 :
1044 44 : subroutine set_Psurf_BC(ierr)
1045 : integer, intent(out) :: ierr
1046 : logical :: test_partials
1047 : type(auto_diff_real_star_order1) :: lnP1_ad
1048 : include 'formats'
1049 : !test_partials = (1 == s% solver_test_partials_k)
1050 44 : test_partials = .false.
1051 44 : ierr = 0
1052 44 : lnP1_ad = wrap_lnPeos_00(s,1)
1053 44 : resid_ad = lnP_bc_ad/lnP1_ad - 1d0
1054 44 : s% equ(i_P_eqn, 1) = resid_ad%val
1055 : if (test_partials) then
1056 : s% solver_test_partials_val = 0
1057 : end if
1058 : call save_eqn_residual_info( &
1059 44 : s, 1, nvar, i_P_eqn, resid_ad, 'set_Psurf_BC', ierr)
1060 : if (test_partials) then
1061 : s% solver_test_partials_var = 0
1062 : s% solver_test_partials_dval_dx = 0
1063 : write(*,*) 'set_Psurf_BC', s% solver_test_partials_var
1064 : end if
1065 44 : end subroutine set_Psurf_BC
1066 :
1067 :
1068 0 : subroutine set_momentum_BC(ierr)
1069 : use hydro_riemann, only: do_surf_Riemann_dudt_eqn
1070 : use hydro_momentum, only: do_surf_momentum_eqn
1071 : integer, intent(out) :: ierr
1072 : include 'formats'
1073 : ierr = 0
1074 0 : if (s% u_flag) then
1075 0 : call do_surf_Riemann_dudt_eqn(s, P_bc_ad, nvar, ierr)
1076 : else
1077 0 : call do_surf_momentum_eqn(s, P_bc_ad, nvar, ierr)
1078 : end if
1079 0 : end subroutine set_momentum_BC
1080 :
1081 :
1082 0 : subroutine set_compression_BC(ierr)
1083 : integer, intent(out) :: ierr
1084 : type(auto_diff_real_star_order1) :: &
1085 : rho1, rho2, dlnd1, dlnd2
1086 : include 'formats'
1087 : ! gradient of compression vanishes fixes density for cell 1
1088 : ! d_dt(1/rho(1)) = d_dt(1/rho(2)) e.g., Grott, Chernigovski, Glatzel, 2005.
1089 0 : ierr = 0
1090 0 : rho1 = wrap_d_00(s,1)
1091 0 : rho2 = wrap_d_p1(s,1)
1092 0 : dlnd1 = wrap_dxh_lnd(s,1) ! lnd(1) - lnd_start(1)
1093 0 : dlnd2 = shift_p1(wrap_dxh_lnd(s,2)) ! lnd(2) - lnd_start(2)
1094 0 : resid_ad = (rho2*dlnd1 - rho1*dlnd2)/s% dt
1095 0 : s% equ(i_P_eqn, 1) = resid_ad%val
1096 : call save_eqn_residual_info( &
1097 0 : s, 1, nvar, i_P_eqn, resid_ad, 'set_compression_BC', ierr)
1098 0 : end subroutine set_compression_BC
1099 :
1100 :
1101 0 : subroutine set_fixed_vsurf_outer_BC(ierr)
1102 : integer, intent(out) :: ierr
1103 : type(auto_diff_real_star_order1) :: vsurf
1104 : include 'formats'
1105 0 : ierr = 0
1106 0 : if (do_du_dt) then
1107 0 : vsurf = wrap_u_00(s,1)
1108 0 : else if (do_dv_dt) then
1109 0 : vsurf = wrap_v_00(s,1)
1110 : else
1111 0 : ierr = -1
1112 0 : write(*,*) 'set_fixed_vsurf_outer_BC requires u_flag or v_flag true'
1113 0 : return
1114 : end if
1115 0 : resid_ad = (vsurf - s% fixed_vsurf)/s% csound_start(1)
1116 0 : s% equ(i_P_eqn,1) = resid_ad%val
1117 : call save_eqn_residual_info( &
1118 0 : s, 1, nvar, i_P_eqn, resid_ad, 'set_fixed_vsurf_outer_BC', ierr)
1119 : end subroutine set_fixed_vsurf_outer_BC
1120 :
1121 :
1122 : end subroutine PT_eqns_surf
1123 :
1124 :
1125 : integer function equ_extra_profile_columns(id)
1126 : use star_def, only: star_info
1127 : integer, intent(in) :: id
1128 : type (star_info), pointer :: s
1129 : integer :: ierr
1130 : ierr = 0
1131 : call star_ptr(id, s, ierr)
1132 : if (ierr /= 0) return
1133 : equ_extra_profile_columns = s% nvar_hydro
1134 : end function equ_extra_profile_columns
1135 :
1136 :
1137 : subroutine equ_data_for_extra_profile_columns( &
1138 : id, n, nz, names, vals, ierr)
1139 : use star_def, only: maxlen_profile_column_name, star_info
1140 : integer, intent(in) :: id, n, nz
1141 : character (len=maxlen_profile_column_name) :: names(n)
1142 : real(dp) :: vals(nz,n)
1143 : integer, intent(out) :: ierr
1144 : integer :: i, k
1145 : type (star_info), pointer :: s
1146 : include 'formats'
1147 : ierr = 0
1148 : call star_ptr(id, s, ierr)
1149 : if (ierr /= 0) return
1150 : if (s% nvar_hydro /= n) then
1151 : write(*,3) 'nvar_hydro /= n', s% nvar_hydro, n
1152 : call mesa_error(__FILE__,__LINE__,'equ_data_for_extra_profile_columns')
1153 : end if
1154 : do i=1,n
1155 : do k=1,nz
1156 : vals(k,i) = s% equ(i,k)
1157 : end do
1158 : names(i) = s% nameofequ(i)
1159 : end do
1160 : end subroutine equ_data_for_extra_profile_columns
1161 :
1162 : end module hydro_eqns
|
