Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2015-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_riemann
21 :
22 : use star_private_def
23 : use const_def, only: dp, pi
24 : use star_utils, only: em1, e00, ep1
25 : use utils_lib
26 : use auto_diff
27 : use auto_diff_support
28 :
29 : implicit none
30 :
31 : ! Cheng, J, Shu, C-W, and Zeng, Q.,
32 : ! "A Conservative Lagrangian Scheme for Solving
33 : ! Compressible Fluid Flows with Multiple Internal Energy Equations",
34 : ! Commun. Comput. Phys., 12, pp 1307-1328, 2012.
35 :
36 : ! Cheng, J. and Shu, C-W,
37 : ! "Positivity-preserving Lagrangian scheme for multi-material
38 : ! compressible flow", J. Comp. Phys., 257 (2014), 143-168.
39 :
40 : ! Kappeli, R. and Mishra, S.,
41 : ! "Well-balanced schemes for the Euler equations with gravitation",
42 : ! J. Comp. Phys., 259 (2014), 199-219.
43 :
44 : private
45 : public :: do_surf_Riemann_dudt_eqn, do1_Riemann_momentum_eqn, &
46 : do_uface_and_Pface
47 : ! Riemann energy eqn is now part of the standard energy equation
48 : ! Riemann dlnR_dt rqn is now part of the standard radius equation
49 :
50 : contains
51 :
52 0 : subroutine do_surf_Riemann_dudt_eqn(s, P_surf_ad, nvar, ierr)
53 : type (star_info), pointer :: s
54 : type(auto_diff_real_star_order1), intent(in) :: P_surf_ad
55 : integer, intent(in) :: nvar
56 : integer, intent(out) :: ierr
57 0 : call do1_dudt_eqn(s, 1, P_surf_ad, nvar, ierr)
58 0 : end subroutine do_surf_Riemann_dudt_eqn
59 :
60 :
61 0 : subroutine do1_Riemann_momentum_eqn(s, k, nvar, ierr)
62 : type (star_info), pointer :: s
63 : integer, intent(in) :: k
64 : integer, intent(in) :: nvar
65 : integer, intent(out) :: ierr
66 : type(auto_diff_real_star_order1) :: P_surf_ad
67 0 : P_surf_ad = 0
68 0 : call do1_dudt_eqn(s, k, P_surf_ad, nvar, ierr)
69 0 : end subroutine do1_Riemann_momentum_eqn
70 :
71 :
72 0 : subroutine do1_dudt_eqn( &
73 : s, k, P_surf_ad, nvar, ierr)
74 : use accurate_sum_auto_diff_star_order1
75 : use star_utils, only: get_area_info_opt_time_center, save_eqn_residual_info
76 : use tdc_hydro, only: compute_tdc_Uq_dm_cell
77 : type (star_info), pointer :: s
78 : integer, intent(in) :: k
79 : type(auto_diff_real_star_order1), intent(in) :: P_surf_ad ! only for k=1
80 : integer, intent(in) :: nvar
81 : integer, intent(out) :: ierr
82 : integer :: nz, i_du_dt
83 : type(auto_diff_real_star_order1) :: &
84 : flux_in_ad, flux_out_ad, diffusion_source_ad, &
85 : geometry_source_ad, gravity_source_ad, &
86 : area_00, area_p1, inv_R2_00, inv_R2_p1, &
87 : dudt_expected_ad, dudt_actual_ad, resid_ad, &
88 : Uq_cell
89 : type(accurate_auto_diff_real_star_order1) :: sum_ad
90 : real(dp) :: dt, dm, ie_plus_ke, scal, residual
91 : logical :: dbg, do_diffusion, test_partials
92 : real(dp) :: v_drag, drag_factor, drag_fraction
93 :
94 : include 'formats'
95 0 : dbg = .false.
96 :
97 : !test_partials = (k == s% solver_test_partials_k)
98 0 : test_partials = .false.
99 :
100 0 : if (s% use_other_momentum) &
101 0 : call mesa_error(__FILE__,__LINE__,'Riemann dudt does not support use_other_momentum')
102 0 : if (s% use_other_momentum_implicit) &
103 0 : call mesa_error(__FILE__,__LINE__,'Riemann dudt does not support use_other_momentum_implicit')
104 0 : if (s% use_mass_corrections) &
105 0 : call mesa_error(__FILE__,__LINE__,'Riemann dudt does not support use_mass_corrections')
106 :
107 : ierr = 0
108 0 : nz = s% nz
109 0 : i_du_dt = s% i_du_dt
110 0 : dt = s% dt
111 0 : dm = s% dm(k)
112 :
113 0 : call get_area_info_opt_time_center(s, k, area_00, inv_R2_00, ierr)
114 0 : if (ierr /= 0) return
115 0 : if (k < nz) then
116 0 : call get_area_info_opt_time_center(s, k+1, area_p1, inv_R2_p1, ierr)
117 0 : if (ierr /= 0) return
118 0 : area_p1 = shift_p1(area_p1)
119 0 : inv_R2_p1 = shift_p1(inv_R2_p1)
120 : end if
121 :
122 0 : call setup_momentum_flux
123 0 : call setup_geometry_source(ierr); if (ierr /= 0) return
124 0 : call setup_gravity_source
125 0 : call setup_diffusion_source
126 :
127 : ! Add turbulent eddy viscous acceleration Uq for TDC as source
128 0 : Uq_cell = 0d0
129 0 : if (s% MLT_option == 'TDC' .and. s%TDC_alpha_M > 0d0) then
130 0 : Uq_cell = compute_tdc_Uq_dm_cell(s, k, ierr) ! Uq * dm
131 0 : if (ierr /= 0) return
132 : end if
133 :
134 0 : sum_ad = flux_in_ad
135 0 : sum_ad = sum_ad - flux_out_ad
136 0 : sum_ad = sum_ad + geometry_source_ad
137 0 : sum_ad = sum_ad + gravity_source_ad
138 0 : sum_ad = sum_ad + diffusion_source_ad
139 0 : sum_ad = sum_ad + Uq_cell
140 0 : dudt_expected_ad = sum_ad
141 0 : dudt_expected_ad = dudt_expected_ad/dm
142 :
143 : ! implement drag
144 0 : drag_factor = s% v_drag_factor
145 0 : v_drag = s% v_drag
146 0 : if (s% q(k) < s% q_for_v_drag_full_off) then
147 : drag_fraction = 0d0
148 0 : else if (s% q(k) > s% q_for_v_drag_full_on) then
149 : drag_fraction = 1d0
150 : else
151 : drag_fraction = (s% q(k) - s% q_for_v_drag_full_off)&
152 0 : /(s% q_for_v_drag_full_on - s% q_for_v_drag_full_off)
153 : end if
154 0 : drag_factor = drag_factor*drag_fraction
155 :
156 0 : if (drag_factor > 0d0) then
157 0 : if (s% u(k) > v_drag) then
158 0 : dudt_expected_ad = dudt_expected_ad - drag_factor*pow2(s% u(k) - v_drag)/s% r(k)
159 0 : else if (s% u(k) < -v_drag) then
160 0 : dudt_expected_ad = dudt_expected_ad + drag_factor*pow2(s% u(k) + v_drag)/s% r(k)
161 : end if
162 : end if
163 :
164 :
165 : ! make residual units be relative difference in energy
166 0 : ie_plus_ke = s% energy_start(k) + 0.5d0*s% u_start(k)*s% u_start(k)
167 0 : scal = dt*max(abs(s% u_start(k)),s% csound_start(k))/ie_plus_ke
168 0 : if (k == 1) scal = scal*1d-2
169 :
170 0 : dudt_actual_ad = 0d0
171 0 : dudt_actual_ad%val = s% dxh_u(k)/dt
172 0 : dudt_actual_ad%d1Array(i_v_00) = 1d0/dt
173 :
174 0 : resid_ad = scal*(dudt_expected_ad - dudt_actual_ad)
175 0 : residual = resid_ad%val
176 0 : s% equ(i_du_dt, k) = residual
177 :
178 0 : if (is_bad(residual)) then
179 0 : ierr = -1
180 0 : return
181 : !$omp critical (dudt_eqn)
182 : write(*,2) 'residual', k, residual
183 : call mesa_error(__FILE__,__LINE__,'do1_dudt_eqn')
184 : !$omp end critical (dudt_eqn)
185 : end if
186 :
187 0 : call save_eqn_residual_info(s, k, nvar, i_du_dt, resid_ad, 'do1_dudt_eqn', ierr)
188 :
189 : if (test_partials) then
190 : s% solver_test_partials_val = resid_ad% val
191 : end if
192 :
193 : if (test_partials) then
194 : s% solver_test_partials_var = s% i_lnR
195 : s% solver_test_partials_dval_dx = resid_ad% d1Array(i_lnR_00)
196 : !write(*,*) 'do1_dudt_eqn', s% solver_test_partials_var
197 : end if
198 :
199 : contains
200 :
201 0 : subroutine setup_momentum_flux
202 0 : if (k == 1) then
203 0 : flux_out_ad = P_surf_ad*area_00
204 : else
205 0 : flux_out_ad = s% P_face_ad(k)*area_00
206 : end if
207 0 : if (k < nz) then
208 0 : flux_in_ad = shift_p1(s% P_face_ad(k+1))*area_p1
209 : else
210 0 : flux_in_ad = 0d0
211 : end if
212 0 : end subroutine setup_momentum_flux
213 :
214 0 : subroutine setup_geometry_source(ierr)
215 : use star_utils, only: calc_Ptot_ad_tw
216 : integer, intent(out) :: ierr
217 : type(auto_diff_real_star_order1) :: P
218 0 : real(dp), dimension(s% species) :: d_Ptot_dxa
219 : logical, parameter :: skip_Peos = .false., skip_mlt_Pturb = .false.
220 : ierr = 0
221 : ! use same P here as the cell pressure in P_face calculation
222 0 : call calc_Ptot_ad_tw(s, k, skip_Peos, skip_mlt_Pturb, P, d_Ptot_dxa, ierr)
223 0 : if (ierr /= 0) return
224 0 : if (k == nz) then
225 : ! no flux in from left, so only have geometry source on right
226 : ! this matters for cases with R_center > 0.
227 0 : geometry_source_ad = P*area_00
228 : else
229 0 : geometry_source_ad = P*(area_00 - area_p1)
230 : end if
231 : end subroutine setup_geometry_source
232 :
233 0 : subroutine setup_gravity_source
234 : type(auto_diff_real_star_order1) :: G00, Gp1, gsL, gsR
235 : real(dp) :: mR, mL
236 : ! left 1/2 of dm gets gravity force at left face
237 : ! right 1/2 of dm gets gravity force at right face.
238 : ! this form is to match the gravity force equilibrium reconstruction.
239 0 : mR = s% m(k)
240 0 : if (k == nz) then
241 0 : mL = s% M_center
242 : else
243 0 : mL = s% m(k+1)
244 : end if
245 0 : call get_G(s, k, G00)
246 0 : gsR = -G00*mR*0.5d0*dm*inv_R2_00
247 0 : if (k == nz) then
248 0 : gsL = 0d0
249 : else
250 0 : call get_G(s, k+1, Gp1)
251 0 : Gp1 = shift_p1(Gp1)
252 0 : gsL = -Gp1*mL*0.5d0*dm*inv_R2_p1
253 : end if
254 0 : gravity_source_ad = gsL + gsR ! total gravitational force on cell
255 :
256 :
257 0 : end subroutine setup_gravity_source
258 :
259 0 : subroutine setup_diffusion_source
260 : type(auto_diff_real_star_order1) :: u_m1, u_00, u_p1
261 : real(dp) :: sig00, sigp1
262 0 : do_diffusion = s% RTI_flag .and. s% dudt_RTI_diffusion_factor > 0d0
263 0 : if (do_diffusion) then ! add diffusion source term to dudt
264 0 : u_p1 = 0d0 ! sets val and d1Array to 0
265 0 : if (k < nz) then
266 0 : sigp1 = s% dudt_RTI_diffusion_factor*s% sig_RTI(k+1)
267 0 : u_p1%val = s% u(k+1)
268 0 : u_p1%d1Array(i_v_p1) = 1d0
269 : else
270 0 : sigp1 = 0
271 : end if
272 0 : u_m1 = 0d0 ! sets val and d1Array to 0
273 0 : if (k > 1) then
274 0 : sig00 = s% dudt_RTI_diffusion_factor*s% sig_RTI(k)
275 0 : u_m1%val = s% u(k-1)
276 0 : u_m1%d1Array(i_v_m1) = 1d0
277 : else
278 0 : sig00 = 0
279 : end if
280 0 : u_00 = 0d0 ! sets val and d1Array to 0
281 0 : u_00%val = s% u(k)
282 0 : u_00%d1Array(i_v_00) = 1d0
283 0 : diffusion_source_ad = sig00*(u_m1 - u_00) - sigp1*(u_00 - u_p1)
284 : else
285 0 : diffusion_source_ad = 0d0
286 : end if
287 0 : s% dudt_RTI(k) = diffusion_source_ad%val/dm
288 0 : end subroutine setup_diffusion_source
289 :
290 : end subroutine do1_dudt_eqn
291 :
292 :
293 0 : subroutine do_uface_and_Pface(s, ierr)
294 : type (star_info), pointer :: s
295 : integer, intent(out) :: ierr
296 : integer :: k, op_err
297 : include 'formats'
298 0 : ierr = 0
299 0 : !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2)
300 : do k = 1, s% nz
301 : op_err = 0
302 : call do1_uface_and_Pface(s, k, op_err)
303 : if (op_err /= 0) ierr = op_err
304 : end do
305 : !$OMP END PARALLEL DO
306 0 : end subroutine do_uface_and_Pface
307 :
308 :
309 0 : subroutine get_G(s, k, G)
310 : type (star_info), pointer :: s
311 : integer, intent(in) :: k
312 : type(auto_diff_real_star_order1), intent(out) :: G
313 : real(dp) :: cgrav
314 0 : cgrav = s% cgrav(k)
315 0 : G = cgrav
316 0 : if (s% rotation_flag .and. s% use_gravity_rotation_correction) &
317 0 : G = G*s% fp_rot(k)
318 0 : end subroutine get_G
319 :
320 :
321 0 : subroutine do1_uface_and_Pface(s, k, ierr)
322 : use eos_def, only: i_gamma1, i_lnfree_e, i_lnPgas
323 : use star_utils, only: calc_Ptot_ad_tw, get_face_weights
324 : use hydro_rsp2, only: compute_Uq_face
325 : type (star_info), pointer :: s
326 : integer, intent(in) :: k
327 : integer, intent(out) :: ierr
328 : logical :: test_partials
329 :
330 : type(auto_diff_real_star_order1) :: &
331 : r_ad, A_ad, PL_ad, PR_ad, uL_ad, uR_ad, rhoL_ad, rhoR_ad, &
332 : gamma1L_ad, gamma1R_ad, csL_ad, csR_ad, G_ad, dPdm_grav_ad, &
333 : Sl1_ad, Sl2_ad, Sr1_ad, Sr2_ad, numerator_ad, denominator_ad, &
334 : Sl_ad, Sr_ad, Ss_ad, P_face_L_ad, P_face_R_ad, du_ad, Uq_ad
335 0 : real(dp), dimension(s% species) :: d_Ptot_dxa ! skip this
336 : logical, parameter :: skip_Peos = .false., skip_mlt_Pturb = .false.
337 : real(dp) :: delta_m, f
338 :
339 : include 'formats'
340 :
341 0 : ierr = 0
342 0 : test_partials = .false.
343 : !test_partials = (k == s% solver_test_partials_k)
344 :
345 0 : s% RTI_du_diffusion_kick(k) = 0d0
346 0 : s% d_uface_domega(k) = 0
347 :
348 0 : if (k == 1) then
349 0 : s% u_face_ad(k) = wrap_u_00(s,k)
350 0 : s% P_face_ad(k) = wrap_Peos_00(s,k)
351 0 : return
352 : end if
353 :
354 0 : r_ad = wrap_r_00(s,k)
355 0 : A_ad = 4d0*pi*pow2(r_ad)
356 :
357 0 : call calc_Ptot_ad_tw(s, k, skip_Peos, skip_mlt_Pturb, PL_ad, d_Ptot_dxa, ierr)
358 0 : if (ierr /= 0) return
359 0 : call calc_Ptot_ad_tw(s, k-1, skip_Peos, skip_mlt_Pturb, PR_ad, d_Ptot_dxa, ierr)
360 0 : if (ierr /= 0) return
361 0 : PR_ad = shift_m1(PR_ad)
362 :
363 0 : uL_ad = wrap_u_00(s,k)
364 0 : uR_ad = wrap_u_m1(s,k)
365 :
366 0 : rhoL_ad = wrap_d_00(s,k)
367 0 : rhoR_ad = wrap_d_m1(s,k)
368 :
369 0 : gamma1L_ad = wrap_gamma1_00(s,k)
370 0 : gamma1R_ad = wrap_gamma1_m1(s,k)
371 :
372 0 : csL_ad = sqrt(gamma1L_ad*PL_ad/rhoL_ad)
373 0 : csR_ad = sqrt(gamma1R_ad*PR_ad/rhoR_ad)
374 :
375 : ! change PR and PL for gravity
376 0 : call get_G(s, k, G_ad)
377 :
378 0 : dPdm_grav_ad = -G_ad*s% m_grav(k)/(pow2(r_ad)*A_ad) ! cm^-1 s^-2
379 :
380 0 : delta_m = 0.5d0*s% dm(k) ! positive delta_m from left center to edge
381 0 : PL_ad = PL_ad + delta_m*dPdm_grav_ad
382 :
383 0 : delta_m = -0.5d0*s% dm(k-1) ! negative delta_m from right center to edge
384 0 : PR_ad = PR_ad + delta_m*dPdm_grav_ad
385 :
386 : ! acoustic wavespeeds (eqn 2.38)
387 0 : Sl1_ad = uL_ad - csL_ad
388 0 : Sl2_ad = uR_ad - csR_ad
389 :
390 : ! take Sl = min(Sl1, Sl2)
391 0 : if (Sl1_ad%val < Sl2_ad%val) then
392 0 : Sl_ad = Sl1_ad
393 : else
394 0 : Sl_ad = Sl2_ad
395 : end if
396 :
397 0 : Sr1_ad = uR_ad + csR_ad
398 0 : Sr2_ad = uL_ad + csL_ad
399 :
400 : ! take Sr = max(Sr1, Sr2)
401 0 : if (Sr1_ad%val > Sr2_ad%val) then
402 0 : Sr_ad = Sr1_ad
403 : else
404 0 : Sr_ad = Sr2_ad
405 : end if
406 :
407 : ! contact velocity (eqn 2.20)
408 0 : numerator_ad = uR_ad*rhoR_ad*(Sr_ad - uR_ad) + uL_ad*rhoL_ad*(uL_ad - Sl_ad) + (PL_ad - PR_ad)
409 0 : denominator_ad = rhoR_ad*(Sr_ad - uR_ad) + rhoL_ad*(uL_ad - Sl_ad)
410 :
411 0 : if (denominator_ad%val == 0d0 .or. is_bad(denominator_ad%val)) then
412 0 : ierr = -1
413 0 : if (s% report_ierr) then
414 0 : write(*,2) 'u_face denominator bad', k, denominator_ad%val
415 : end if
416 0 : return
417 : end if
418 :
419 0 : Ss_ad = numerator_ad/denominator_ad
420 :
421 0 : s% u_face_ad(k) = Ss_ad
422 0 : s% d_uface_domega(k) = s% u_face_ad(k)%d1Array(i_L_00)
423 :
424 : ! contact pressure (eqn 2.19)
425 0 : P_face_L_ad = rhoL_ad*(uL_ad-Sl_ad)*(uL_ad-Ss_ad) + PL_ad
426 0 : P_face_R_ad = rhoR_ad*(uR_ad-Sr_ad)*(uR_ad-Ss_ad) + PR_ad
427 :
428 0 : s% P_face_ad(k) = 0.5d0*(P_face_L_ad + P_face_R_ad) ! these are ideally equal
429 :
430 0 : if (k < s% nz .and. s% RTI_flag) then
431 : if (s% eta_RTI(k) > 0d0 .and. &
432 0 : s% dlnddt_RTI_diffusion_factor > 0d0 .and. s% dt > 0d0) then
433 0 : f = s% dlnddt_RTI_diffusion_factor*s% eta_RTI(k)/s% dm_bar(k)
434 0 : du_ad = f*A_ad*(rhoL_ad - rhoR_ad) ! bump uface in direction of lower density
435 0 : s% RTI_du_diffusion_kick(k) = du_ad%val
436 0 : s% u_face_ad(k) = s% u_face_ad(k) + du_ad
437 : end if
438 : end if
439 :
440 :
441 0 : if (s% RSP2_flag) then ! include Uq in u_face, To do: implement in sources instead ~ EbF
442 0 : Uq_ad = compute_Uq_face(s, k, ierr)
443 0 : if (ierr /= 0) return
444 0 : s% u_face_ad(k) = s% u_face_ad(k) + Uq_ad
445 : end if
446 :
447 0 : s% u_face_val(k) = s% u_face_ad(k)%val
448 :
449 0 : if (s% P_face_start(k) < 0d0) then
450 0 : s% u_face_start(k) = s% u_face_val(k)
451 0 : s% P_face_start(k) = s% P_face_ad(k)%val
452 : end if
453 :
454 : if (test_partials) then
455 : s% solver_test_partials_val = PL_ad% val
456 : s% solver_test_partials_var = s% i_w_div_wc
457 : s% solver_test_partials_dval_dx = PL_ad% d1Array(i_w_div_wc_00)
458 : write(*,*) 'do1_uface_and_Pface', s% solver_test_partials_var, PL_ad% val
459 : end if
460 :
461 : end subroutine do1_uface_and_Pface
462 :
463 : end module hydro_riemann
|
