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_momentum
21 :
22 : use star_private_def
23 : use const_def, only: dp, ln10, secyer
24 : use utils_lib, only: mesa_error, is_bad
25 : use auto_diff
26 : use star_utils, only: em1, e00, ep1
27 :
28 : implicit none
29 :
30 : private
31 : public :: do1_momentum_eqn
32 : public :: do_surf_momentum_eqn
33 : public :: do1_radius_eqn
34 : public :: expected_HSE_grav_term
35 :
36 : contains
37 :
38 0 : subroutine do_surf_momentum_eqn(s, P_surf_ad, nvar, ierr)
39 : use star_utils, only: store_partials
40 : type (star_info), pointer :: s
41 : type(auto_diff_real_star_order1), intent(in) :: P_surf_ad
42 : integer, intent(in) :: nvar
43 : integer, intent(out) :: ierr
44 0 : real(dp) :: d_dm1(nvar), d_d00(nvar), d_dp1(nvar)
45 : include 'formats'
46 : ierr = 0
47 : call get1_momentum_eqn( &
48 0 : s, 1, P_surf_ad, nvar, d_dm1, d_d00, d_dp1, ierr)
49 0 : if (ierr /= 0) then
50 0 : if (s% report_ierr) write(*,2) 'ierr /= 0 for do_surf_momentum_eqn'
51 : return
52 : end if
53 : call store_partials( &
54 0 : s, 1, s% i_dv_dt, nvar, d_dm1, d_d00, d_dp1, 'do_surf_momentum_eqn', ierr)
55 0 : end subroutine do_surf_momentum_eqn
56 :
57 :
58 52304 : subroutine do1_momentum_eqn(s, k, nvar, ierr)
59 : use star_utils, only: store_partials
60 : type (star_info), pointer :: s
61 : integer, intent(in) :: k
62 : integer, intent(in) :: nvar
63 : integer, intent(out) :: ierr
64 52304 : real(dp) :: d_dm1(nvar), d_d00(nvar), d_dp1(nvar)
65 : type(auto_diff_real_star_order1) :: P_surf_ad ! only used if k == 1
66 : include 'formats'
67 52304 : P_surf_ad = 0d0
68 : call get1_momentum_eqn( &
69 52304 : s, k, P_surf_ad, nvar, d_dm1, d_d00, d_dp1, ierr)
70 52304 : if (ierr /= 0) then
71 0 : if (s% report_ierr) write(*,2) 'ierr /= 0 for get1_momentum_eqn', k
72 : return
73 : end if
74 : call store_partials( &
75 52304 : s, k, s% i_dv_dt, nvar, d_dm1, d_d00, d_dp1, 'do1_momentum_eqn', ierr)
76 : end subroutine do1_momentum_eqn
77 :
78 :
79 52304 : subroutine get1_momentum_eqn( &
80 : s, k, P_surf_ad, nvar, &
81 52304 : d_dm1, d_d00, d_dp1, ierr)
82 : use accurate_sum_auto_diff_star_order1
83 : use auto_diff_support
84 :
85 : type (star_info), pointer :: s
86 : integer, intent(in) :: k
87 : type(auto_diff_real_star_order1), intent(in) :: P_surf_ad ! only used if k == 1
88 : integer, intent(in) :: nvar
89 : real(dp), intent(out) :: d_dm1(nvar), d_d00(nvar), d_dp1(nvar)
90 : integer, intent(out) :: ierr
91 :
92 : real(dp) :: residual, dm_face, dPtot, iPtotavg, dm_div_A
93 : real(dp), dimension(s% species) :: &
94 104608 : d_dPtot_dxam1, d_dPtot_dxa00, d_iPtotavg_dxam1, d_iPtotavg_dxa00, &
95 52304 : d_residual_dxam1, d_residual_dxa00
96 : integer :: nz, i_dv_dt, i_lum, i_v
97 : logical :: test_partials
98 :
99 : type(auto_diff_real_star_order1) :: resid1_ad, resid_ad, &
100 : other_ad, dm_div_A_ad, grav_ad, area_ad, dPtot_ad, d_mlt_Pturb_ad, &
101 : iPtotavg_ad, other_dm_div_A_ad, grav_dm_div_A_ad, &
102 : RTI_terms_ad, RTI_terms_dm_div_A_ad, accel_ad, Uq_ad
103 : type(accurate_auto_diff_real_star_order1) :: residual_sum_ad
104 :
105 : include 'formats'
106 :
107 : !test_partials = (k == s% solver_test_partials_k)
108 52304 : test_partials = .false.
109 :
110 52304 : ierr = 0
111 52304 : call init
112 :
113 : ! dv/dt = - G*m/r^2 - (dPtot_ad + d_mlt_Pturb_ad)*area/dm + extra_grav + Uq + RTI_diffusion + RTI_kick
114 : !
115 : ! grav_ad = expected_HSE_grav_term = -G*m/r^2 with possible modifications for rotation
116 : ! other_ad = expected_non_HSE_term = extra_grav - dv/dt + Uq
117 : ! extra_grav is from the other_momentum hook
118 : ! dPtot_ad = pressure difference across face from center to center of adjacent cells (excluding mlt_Pturb effects)
119 : ! Ptot = P_ad + Pvsc_ad + Ptrb_ad + extra_pressure, with time centering
120 : ! iPtotavg_ad = 1/(avg Ptot). for normalizing equation
121 : ! d_mlt_Pturb_ad = difference in MLT convective pressure across face
122 : ! RTI_terms_ad = RTI_diffusion + RTI_kick
123 : ! dm_div_A_ad = dm/area
124 : !
125 : ! 0 = extra_grav - dv/dt + Uq - G*m/r^2 - RTI_diffusion - RTI_kick - (dPtot_ad + d_mlt_Pturb_ad)*area/dm
126 : ! 0 = other + grav - RTI_terms - (dPtot_ad + d_mlt_Pturb_ad)*area/dm
127 : ! 0 = (other + grav - RTI_terms)*dm/area - dPtot_ad - d_mlt_Pturb_ad
128 : ! 0 = other_dm_div_A_ad + grav_dm_div_A_ad - dPtot_ad - d_mlt_Pturb_ad + RTI_terms_dm_div_A_ad
129 :
130 52304 : call setup_HSE(dm_div_A, ierr); if (ierr /= 0) return ! grav_ad and dm_div_A_ad
131 52304 : call setup_non_HSE(ierr); if (ierr /= 0) return ! other = s% extra_grav(k) - dv_dt
132 52304 : call setup_dPtot(ierr); if (ierr /= 0) return ! dPtot_ad, iPtotavg_ad
133 52304 : call setup_d_mlt_Pturb(ierr); if (ierr /= 0) return ! d_mlt_Pturb_ad
134 52304 : call setup_RTI_terms(ierr); if (ierr /= 0) return ! RTI_terms_ad
135 :
136 52304 : other_dm_div_A_ad = other_ad*dm_div_A_ad
137 52304 : grav_dm_div_A_ad = grav_ad*dm_div_A_ad
138 52304 : RTI_terms_dm_div_A_ad = RTI_terms_ad*dm_div_A_ad
139 :
140 : ! sum terms in residual_sum_ad using accurate_auto_diff_real_star_order1
141 52304 : residual_sum_ad = other_dm_div_A_ad
142 52304 : residual_sum_ad = residual_sum_ad + grav_dm_div_A_ad
143 52304 : residual_sum_ad = residual_sum_ad - dPtot_ad
144 52304 : residual_sum_ad = residual_sum_ad - d_mlt_Pturb_ad
145 52304 : residual_sum_ad = residual_sum_ad + RTI_terms_dm_div_A_ad
146 :
147 : resid1_ad = residual_sum_ad ! convert back to auto_diff_real_star_order1
148 52304 : resid_ad = resid1_ad*iPtotavg_ad ! scaling
149 52304 : residual = resid_ad%val
150 52304 : s% equ(i_dv_dt, k) = residual
151 :
152 : !s% xtra1_array(k) = s% Peos(k)
153 : !s% xtra2_array(k) = 1d0/s% rho(k)
154 : !s% xtra3_array(k) = s% T(k)
155 : !s% xtra4_array(k) = s% v(k)
156 : !s% xtra5_array(k) = s% etrb(k)
157 : !s% xtra6_array(k) = s% r(k)
158 :
159 52304 : if (is_bad(residual)) then
160 0 : !$omp critical (hydro_momentum_crit1)
161 0 : write(*,2) 'momentum eqn residual', k, residual
162 0 : call mesa_error(__FILE__,__LINE__,'get1_momentum_eqn')
163 : !$omp end critical (hydro_momentum_crit1)
164 : end if
165 : if (test_partials) then
166 : s% solver_test_partials_val = residual
167 : end if
168 52304 : call unpack_res18(s% species, resid_ad)
169 :
170 : if (test_partials) then
171 : s% solver_test_partials_var = 0
172 : s% solver_test_partials_dval_dx = d_d00(s% solver_test_partials_var)
173 : write(*,*) 'get1_momentum_eqn', s% solver_test_partials_var
174 : end if
175 :
176 : contains
177 :
178 52304 : subroutine init
179 52304 : i_dv_dt = s% i_dv_dt
180 52304 : i_lum = s% i_lum
181 52304 : i_v = s% i_v
182 52304 : nz = s% nz
183 : ! in the momentum equation, e.g., dP/dr = -g * rho (for HSE),
184 : ! rho represents the inertial (gravitational) mass density.
185 : ! since dm is baryonic mass, correct dm_face when using mass corrections
186 : ! this will be used in the calculation of dm_div_A
187 52304 : if (s% use_mass_corrections) then
188 0 : if (k > 1) then
189 0 : dm_face = (s% dm(k)*s% mass_correction(k) + s% dm(k-1)*s% mass_correction(k-1))/2d0
190 : else ! k == 1
191 0 : dm_face = s% dm(k)*s% mass_correction(k)/2d0
192 : end if
193 : else
194 52304 : if (k > 1) then
195 52304 : dm_face = (s% dm(k) + s% dm(k-1))/2d0
196 : else ! k == 1
197 0 : dm_face = s% dm(k)/2d0
198 : end if
199 : end if
200 1935248 : d_dm1 = 0d0; d_d00 = 0d0; d_dp1 = 0d0
201 52304 : end subroutine init
202 :
203 52304 : subroutine setup_HSE(dm_div_A, ierr)
204 : real(dp), intent(out) :: dm_div_A
205 : integer, intent(out) :: ierr
206 : include 'formats'
207 : ierr = 0
208 52304 : call expected_HSE_grav_term(s, k, grav_ad, area_ad, ierr)
209 52304 : if (ierr /= 0) return
210 52304 : dm_div_A_ad = dm_face/area_ad
211 52304 : dm_div_A = dm_div_A_ad%val
212 : end subroutine setup_HSE
213 :
214 52304 : subroutine setup_non_HSE(ierr)
215 : integer, intent(out) :: ierr
216 : real(dp) :: other
217 : include 'formats'
218 : ierr = 0
219 : ! other = extra_grav - dv/dt
220 52304 : call expected_non_HSE_term(s, k, other_ad, other, accel_ad, Uq_ad, ierr)
221 52304 : end subroutine setup_non_HSE
222 :
223 52304 : subroutine setup_dPtot(ierr)
224 : integer, intent(out) :: ierr
225 : include 'formats'
226 : ierr = 0
227 : ! dPtot = pressure difference across face from center to center of adjacent cells.
228 : ! iPtotavg = average pressure at face for normalization of the equation to something like dlnP/dm
229 : call get_dPtot_face_info(s, k, P_surf_ad, &
230 : dPtot_ad, dPtot, d_dPtot_dxam1, d_dPtot_dxa00, &
231 52304 : iPtotavg_ad, iPtotavg, d_iPtotavg_dxam1, d_iPtotavg_dxa00, ierr)
232 : if (ierr /= 0) return
233 : end subroutine setup_dPtot
234 :
235 52304 : subroutine setup_d_mlt_Pturb(ierr)
236 : use star_utils, only: get_rho_face
237 : integer, intent(out) :: ierr
238 : type(auto_diff_real_star_order1) :: rho_00, rho_m1
239 52304 : ierr = 0
240 : ! d_mlt_Pturb = difference in MLT convective pressure across face
241 52304 : if (s% mlt_Pturb_factor > 0d0 .and. s% mlt_vc_old(k) > 0d0) then
242 0 : rho_00 = wrap_d_00(s,k)
243 0 : rho_m1 = wrap_d_m1(s,k)
244 0 : d_mlt_Pturb_ad = s% mlt_Pturb_factor*pow2(s% mlt_vc_old(k))*(rho_m1 - rho_00)/3d0
245 : else
246 52304 : d_mlt_Pturb_ad = 0d0
247 : end if
248 52304 : end subroutine setup_d_mlt_Pturb
249 :
250 52304 : subroutine setup_RTI_terms(ierr)
251 : use auto_diff_support
252 : integer, intent(out) :: ierr
253 : type(auto_diff_real_star_order1) :: v_p1, v_00, v_m1, dvdt_diffusion, &
254 : f, rho_00, rho_m1, dvdt_kick
255 : real(dp) :: sigm1, sig00
256 52304 : ierr = 0
257 52304 : RTI_terms_ad = 0d0
258 52304 : if (.not. s% RTI_flag) return
259 0 : if (k >= s% nz .or. k <= 1) return
260 : ! diffusion of specific momentum (i.e. v)
261 0 : if (s% dudt_RTI_diffusion_factor > 0d0) then ! add diffusion source term to dvdt
262 : ! sigmid_RTI(k) is mixing flow at center k in (gm sec^1)
263 0 : sigm1 = s% dudt_RTI_diffusion_factor*s% sigmid_RTI(k-1)
264 0 : sig00 = s% dudt_RTI_diffusion_factor*s% sigmid_RTI(k)
265 0 : v_p1 = wrap_v_p1(s, k)
266 0 : v_00 = wrap_v_00(s, k)
267 0 : v_m1 = wrap_v_m1(s, k)
268 0 : dvdt_diffusion = sig00*(v_p1 - v_00) - sigm1*(v_00 - v_m1) ! (g/s)*(cm/s)
269 0 : dvdt_diffusion = dvdt_diffusion/s% dm_bar(k) ! divide by g to get units of cm/s^2
270 : else
271 0 : dvdt_diffusion = 0d0
272 : end if
273 : ! kick to adjust densities
274 : if (s% eta_RTI(k) > 0d0 .and. &
275 0 : s% dlnddt_RTI_diffusion_factor > 0d0 .and. s% dt > 0d0) then
276 0 : f = s% dlnddt_RTI_diffusion_factor*s% eta_RTI(k)/dm_div_A_ad
277 0 : rho_00 = wrap_d_00(s, k)
278 0 : rho_m1 = wrap_d_m1(s, k)
279 0 : dvdt_kick = f*(rho_00 - rho_m1)/s% dt ! change v according to direction of lower density
280 : else
281 0 : dvdt_kick = 0d0
282 : end if
283 0 : RTI_terms_ad = dvdt_diffusion + dvdt_kick
284 : end subroutine setup_RTI_terms
285 :
286 52304 : subroutine unpack_res18(species, res18)
287 : use star_utils, only: save_eqn_dxa_partials, unpack_residual_partials
288 : integer, intent(in) :: species
289 : type(auto_diff_real_star_order1) :: res18
290 104608 : real(dp) :: resid1, dxap1(species)
291 : logical, parameter :: checking = .true.
292 : integer :: j
293 : include 'formats'
294 : ! do partials wrt composition
295 52304 : resid1 = resid1_ad%val
296 470736 : do j=1,species
297 418432 : d_residual_dxa00(j) = resid1*d_iPtotavg_dxa00(j) - iPtotavg*d_dPtot_dxa00(j)
298 418432 : if (checking) call check_dequ(d_dPtot_dxa00(j),'d_dPtot_dxa00(j)')
299 470736 : if (checking) call check_dequ(d_iPtotavg_dxa00(j),'d_iPtotavg_dxa00(j)')
300 : end do
301 52304 : if (k > 1) then
302 470736 : do j=1,species
303 418432 : d_residual_dxam1(j) = resid1*d_iPtotavg_dxam1(j) - iPtotavg*d_dPtot_dxam1(j)
304 418432 : if (checking) call check_dequ(d_dPtot_dxam1(j),'d_dPtot_dxam1(j)')
305 470736 : if (checking) call check_dequ(d_iPtotavg_dxam1(j),'d_iPtotavg_dxam1(j)')
306 : end do
307 : else
308 0 : d_residual_dxam1 = 0d0
309 : end if
310 470736 : dxap1 = 0d0
311 : call save_eqn_dxa_partials(&
312 : s, k, nvar, i_dv_dt, species, &
313 52304 : d_residual_dxam1, d_residual_dxa00, dxap1, 'get1_momentum_eqn', ierr)
314 : call unpack_residual_partials(s, k, nvar, i_dv_dt, &
315 52304 : res18, d_dm1, d_d00, d_dp1)
316 52304 : end subroutine unpack_res18
317 :
318 1673728 : subroutine check_dequ(dequ, str)
319 : real(dp), intent(in) :: dequ
320 : character (len=*), intent(in) :: str
321 : include 'formats'
322 1673728 : if (is_bad(dequ)) then
323 0 : !$omp critical (hydro_momentum_crit2)
324 0 : ierr = -1
325 0 : if (s% report_ierr) then
326 0 : write(*,2) 'get1_momentum_eqn: bad ' // trim(str), k, dequ
327 : end if
328 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'get1_momentum_eqn')
329 : !$omp end critical (hydro_momentum_crit2)
330 0 : return
331 : end if
332 : end subroutine check_dequ
333 :
334 : end subroutine get1_momentum_eqn
335 :
336 :
337 : ! returns -G*m/r^2 with possible modifications for rotation. MESA 2, eqn 22.
338 104608 : subroutine expected_HSE_grav_term(s, k, grav, area, ierr)
339 : use star_utils, only: get_area_info_opt_time_center
340 : type (star_info), pointer :: s
341 : integer, intent(in) :: k
342 : type(auto_diff_real_star_order1), intent(out) :: area, grav
343 : integer, intent(out) :: ierr
344 :
345 : type(auto_diff_real_star_order1) :: inv_R2
346 : logical :: test_partials
347 :
348 : include 'formats'
349 : ierr = 0
350 :
351 104608 : call get_area_info_opt_time_center(s, k, area, inv_R2, ierr)
352 104608 : if (ierr /= 0) return
353 :
354 104608 : if (s% rotation_flag .and. s% use_gravity_rotation_correction) then
355 0 : grav = -s% cgrav(k)*s% m_grav(k)*inv_R2*s% fp_rot(k)
356 : else
357 104608 : grav = -s% cgrav(k)*s% m_grav(k)*inv_R2
358 : end if
359 :
360 : !test_partials = (k == s% solver_test_partials_k)
361 104608 : test_partials = .false.
362 :
363 : if (test_partials) then
364 : s% solver_test_partials_val = 0
365 : s% solver_test_partials_var = 0
366 : s% solver_test_partials_dval_dx = 0
367 : write(*,*) 'expected_HSE_grav_term', s% solver_test_partials_var
368 : end if
369 :
370 : end subroutine expected_HSE_grav_term
371 :
372 :
373 : ! other = s% extra_grav(k) - s% dv_dt(k)
374 156912 : subroutine expected_non_HSE_term( &
375 : s, k, other_ad, other, accel_ad, Uq_ad, ierr)
376 : use hydro_rsp2, only: compute_Uq_face
377 : use tdc_hydro, only: compute_tdc_Uq_face
378 : use accurate_sum_auto_diff_star_order1
379 : use auto_diff_support
380 : type (star_info), pointer :: s
381 : integer, intent(in) :: k
382 : type(auto_diff_real_star_order1), intent(out) :: &
383 : other_ad, accel_ad,Uq_ad
384 : real(dp), intent(out) :: other
385 : integer, intent(out) :: ierr
386 : type(auto_diff_real_star_order1) :: extra_ad, v_00, &
387 : drag
388 : type(accurate_auto_diff_real_star_order1) :: other_sum_ad
389 : real(dp) :: accel, d_accel_dv
390 : logical :: test_partials, local_v_flag
391 :
392 : include 'formats'
393 :
394 52304 : ierr = 0
395 :
396 52304 : extra_ad = 0d0
397 52304 : if (s% use_other_momentum .or. s% use_other_momentum_implicit) then
398 0 : extra_ad = s% extra_grav(k)
399 : end if
400 :
401 52304 : Uq_ad = 0d0
402 52304 : accel_ad = 0d0
403 52304 : drag = 0d0
404 52304 : s% dvdt_drag(k) = 0d0
405 52304 : if (s% v_flag) then
406 :
407 0 : if (s% i_lnT == 0) then
408 : local_v_flag = .true.
409 : else
410 : local_v_flag = &
411 0 : (s% xh_old(s% i_lnT,k)/ln10 >= s% velocity_logT_lower_bound)
412 : end if
413 :
414 0 : if (local_v_flag) then
415 0 : accel = s% dxh_v(k)/s% dt
416 0 : d_accel_dv = 1d0/s% dt
417 : else ! assume vstart(k) = 0 and
418 : ! constant acceleration dv_dt so vfinal(k) = dv_dt*dt
419 : ! v(k) = dr/dt = average velocity =
420 : ! (vstart + vfinal)/2 = dv_dt*dt/2 when vstart = 0
421 : ! so (1/2)*dv_dt*dt = v(k)
422 0 : accel = 2d0*s% v(k)/s% dt
423 0 : d_accel_dv = 2d0/s% dt
424 : end if
425 0 : accel_ad%val = accel
426 0 : accel_ad%d1Array(i_v_00) = d_accel_dv
427 :
428 0 : if (s% q(k) > s% min_q_for_drag .and. s% drag_coefficient > 0) then
429 0 : v_00 = wrap_v_00(s,k)
430 0 : drag = -s% drag_coefficient*v_00/s% dt
431 0 : s% dvdt_drag(k) = drag%val
432 : end if
433 :
434 0 : if (s% RSP2_flag) then ! Uq(k) is turbulent viscosity drag at face k
435 0 : Uq_ad = compute_Uq_face(s, k, ierr)
436 0 : if (ierr /= 0) return
437 0 : else if (s% TDC_alpha_M > 0 .and. s% MLT_option == 'TDC') then ! Uq(k) is turbulent viscosity drag at face k
438 0 : Uq_ad = compute_tdc_Uq_face(s, k, ierr)
439 0 : if (ierr /= 0) return
440 : end if
441 :
442 : end if ! v_flag
443 :
444 52304 : other_sum_ad = extra_ad
445 52304 : other_sum_ad = other_sum_ad - accel_ad
446 52304 : other_sum_ad = other_sum_ad + drag
447 52304 : other_sum_ad = other_sum_ad + Uq_ad
448 52304 : other_ad = other_sum_ad
449 52304 : other = other_ad%val
450 :
451 : !test_partials = (k == s% solver_test_partials_k)
452 52304 : test_partials = .false.
453 :
454 : if (test_partials) then
455 : s% solver_test_partials_val = 0
456 : s% solver_test_partials_var = 0
457 : s% solver_test_partials_dval_dx = 0d0
458 : write(*,*) 'expected_non_HSE_term', s% solver_test_partials_var
459 : end if
460 :
461 : end subroutine expected_non_HSE_term
462 :
463 : ! dPtot = pressure difference across face from center to center of adjacent cells.
464 : ! excluding mlt_Pturb effects
465 52304 : subroutine get_dPtot_face_info(s, k, P_surf_ad, &
466 52304 : dPtot_ad, dPtot, d_dPtot_dxam1, d_dPtot_dxa00, &
467 52304 : iPtotavg_ad, iPtotavg, d_iPtotavg_dxam1, d_iPtotavg_dxa00, ierr)
468 : use star_utils, only: calc_Ptot_ad_tw
469 : use auto_diff_support
470 : type (star_info), pointer :: s
471 : integer, intent(in) :: k
472 : type(auto_diff_real_star_order1), intent(in) :: P_surf_ad ! only used if k == 1
473 : type(auto_diff_real_star_order1), intent(out) :: dPtot_ad, iPtotavg_ad
474 : real(dp), intent(out) :: dPtot, iPtotavg
475 : real(dp), intent(out), dimension(s% species) :: &
476 : d_dPtot_dxam1, d_dPtot_dxa00, d_iPtotavg_dxam1, d_iPtotavg_dxa00
477 : integer, intent(out) :: ierr
478 :
479 : real(dp) :: Ptotm1, Ptot00, Ptotavg, alfa, beta
480 : real(dp), dimension(s% species) :: &
481 0 : d_Ptotm1_dxam1, d_Ptot00_dxa00, d_Ptotavg_dxam1, d_Ptotavg_dxa00
482 : type(auto_diff_real_star_order1) :: &
483 : Ptot00_ad, Ptotm1_ad, Ptotavg_ad
484 : integer :: j
485 : logical, parameter :: skip_P = .false., skip_mlt_Pturb = .true.
486 : logical :: test_partials
487 :
488 : include 'formats'
489 :
490 : ierr = 0
491 :
492 : call calc_Ptot_ad_tw( &
493 52304 : s, k, skip_P, skip_mlt_Pturb, Ptot00_ad, d_Ptot00_dxa00, ierr)
494 52304 : if (ierr /= 0) return
495 52304 : Ptot00 = Ptot00_ad%val
496 :
497 52304 : if (k > 1) then
498 : call calc_Ptot_ad_tw( &
499 52304 : s, k-1, skip_P, skip_mlt_Pturb, Ptotm1_ad, d_Ptotm1_dxam1, ierr)
500 52304 : if (ierr /= 0) return
501 52304 : Ptotm1_ad = shift_m1(Ptotm1_ad)
502 : else ! k == 1
503 0 : Ptotm1_ad = P_surf_ad
504 : end if
505 52304 : Ptotm1 = Ptotm1_ad%val
506 :
507 52304 : dPtot_ad = Ptotm1_ad - Ptot00_ad
508 52304 : dPtot = Ptotm1 - Ptot00
509 :
510 470736 : do j=1,s% species
511 418432 : d_dPtot_dxam1(j) = d_Ptotm1_dxam1(j)
512 470736 : d_dPtot_dxa00(j) = -d_Ptot00_dxa00(j)
513 : end do
514 :
515 52304 : if (k == 1) then
516 0 : Ptotavg_ad = Ptot00_ad
517 0 : do j=1,s% species
518 0 : d_Ptotavg_dxam1(j) = 0d0
519 0 : d_Ptotavg_dxa00(j) = d_Ptot00_dxa00(j)
520 : end do
521 : else
522 52304 : alfa = s% dq(k-1)/(s% dq(k-1) + s% dq(k))
523 52304 : beta = 1d0 - alfa
524 52304 : Ptotavg_ad = alfa*Ptot00_ad + beta*Ptotm1_ad
525 470736 : do j=1,s% species
526 418432 : d_Ptotavg_dxam1(j) = beta*d_Ptotm1_dxam1(j)
527 470736 : d_Ptotavg_dxa00(j) = alfa*d_Ptot00_dxa00(j)
528 : end do
529 : end if
530 52304 : Ptotavg = Ptotavg_ad%val
531 :
532 52304 : iPtotavg_ad = 1d0/Ptotavg_ad
533 52304 : iPtotavg = 1d0/Ptotavg
534 470736 : do j=1,s% species
535 418432 : d_iPtotavg_dxam1(j) = -iPtotavg*d_Ptotavg_dxam1(j)/Ptotavg
536 470736 : d_iPtotavg_dxa00(j) = -iPtotavg*d_Ptotavg_dxa00(j)/Ptotavg
537 : end do
538 :
539 : !test_partials = (k == s% solver_test_partials_k)
540 52304 : test_partials = .false.
541 :
542 : if (test_partials) then
543 : s% solver_test_partials_val = Ptot00
544 : s% solver_test_partials_var = s% i_lnT
545 : s% solver_test_partials_dval_dx = 0d0
546 : write(*,*) 'get_dPtot_face_info', s% solver_test_partials_var
547 : end if
548 :
549 : end subroutine get_dPtot_face_info
550 :
551 :
552 0 : subroutine do1_radius_eqn(s, k, nvar, ierr)
553 : use auto_diff_support
554 : use star_utils, only: save_eqn_residual_info
555 : type (star_info), pointer :: s
556 : integer, intent(in) :: k, nvar
557 : integer, intent(out) :: ierr
558 : type(auto_diff_real_star_order1) :: &
559 : v00, dxh_lnR, resid_ad, &
560 : dr_div_r0_actual, dr_div_r0_expected
561 : logical :: test_partials, force_zero_v
562 : include 'formats'
563 : !test_partials = (k == s% solver_test_partials_k)
564 0 : test_partials = .false.
565 0 : ierr = 0
566 0 : if (.not. (s% u_flag .or. s% v_flag)) call mesa_error(__FILE__,__LINE__,'must have either v or u for do1_radius_eqn')
567 :
568 : force_zero_v = (s% q(k) > s% velocity_q_upper_bound) .or. &
569 : (s% tau(k) < s% velocity_tau_lower_bound) .or. &
570 : (s% lnT_start(k)/ln10 < s% velocity_logT_lower_bound .and. &
571 0 : s% dt < secyer*s% max_dt_yrs_for_velocity_logT_lower_bound)
572 0 : if (force_zero_v) then
573 0 : if (s% u_flag) then
574 0 : v00 = wrap_u_00(s,k)
575 : else
576 0 : v00 = wrap_v_00(s,k)
577 : end if
578 0 : resid_ad = v00/s% csound_start(k)
579 : call save_eqn_residual_info( &
580 0 : s, k, nvar, s% i_dlnR_dt, resid_ad, 'do1_radius_eqn', ierr)
581 0 : return
582 : end if
583 :
584 : ! dr = r - r0 = v00*dt
585 : ! eqn: dr/r0 = v00*dt/r0
586 : ! (r - r0)/r0 = r/r0 - 1 = exp(lnR)/exp(lnR0) - 1
587 : ! = exp(lnR - lnR0) - 1 = exp(dlnR) - 1 = exp(dlnR_dt*dt) - 1
588 : ! eqn becomes: v00*dt/r0 = expm1(dlnR)
589 0 : dxh_lnR = wrap_dxh_lnR(s,k) ! lnR - lnR_start
590 0 : dr_div_r0_actual = expm1(dxh_lnR) ! expm1(x) = E^x - 1
591 :
592 0 : v00 = wrap_opt_time_center_v_00(s,k)
593 0 : dr_div_r0_expected = v00*s% dt/s% r_start(k)
594 0 : resid_ad = dr_div_r0_expected - dr_div_r0_actual
595 :
596 0 : s% equ(s% i_dlnR_dt, k) = resid_ad%val
597 :
598 : if (test_partials) then
599 : s% solver_test_partials_val = 0
600 : end if
601 : call save_eqn_residual_info( &
602 0 : s, k, nvar, s% i_dlnR_dt, resid_ad, 'do1_radius_eqn', ierr)
603 : if (test_partials) then
604 : s% solver_test_partials_var = 0
605 : s% solver_test_partials_dval_dx = 0
606 : write(*,*) 'do1_radius_eqn', s% solver_test_partials_var
607 : end if
608 : end subroutine do1_radius_eqn
609 :
610 : end module hydro_momentum
|
