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