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 0 : real(dp) :: dt, dm, ie_plus_ke, scal, residual
91 : logical :: dbg, do_diffusion, test_partials
92 0 : 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 : sum_ad = flux_in_ad - flux_out_ad + &
135 0 : geometry_source_ad + gravity_source_ad + diffusion_source_ad + Uq_cell
136 0 : dudt_expected_ad = sum_ad
137 0 : dudt_expected_ad = dudt_expected_ad/dm
138 :
139 : ! implement drag
140 0 : drag_factor = s% v_drag_factor
141 0 : v_drag = s% v_drag
142 0 : if (s% q(k) < s% q_for_v_drag_full_off) then
143 : drag_fraction = 0d0
144 0 : else if (s% q(k) > s% q_for_v_drag_full_on) then
145 : drag_fraction = 1d0
146 : else
147 : drag_fraction = (s% q(k) - s% q_for_v_drag_full_off)&
148 0 : /(s% q_for_v_drag_full_on - s% q_for_v_drag_full_off)
149 : end if
150 0 : drag_factor = drag_factor*drag_fraction
151 :
152 0 : if (drag_factor > 0d0) then
153 0 : if (s% u(k) > v_drag) then
154 0 : dudt_expected_ad = dudt_expected_ad - drag_factor*pow2(s% u(k) - v_drag)/s% r(k)
155 0 : else if (s% u(k) < -v_drag) then
156 0 : dudt_expected_ad = dudt_expected_ad + drag_factor*pow2(s% u(k) + v_drag)/s% r(k)
157 : end if
158 : end if
159 :
160 :
161 : ! make residual units be relative difference in energy
162 0 : ie_plus_ke = s% energy_start(k) + 0.5d0*s% u_start(k)*s% u_start(k)
163 0 : scal = dt*max(abs(s% u_start(k)),s% csound_start(k))/ie_plus_ke
164 0 : if (k == 1) scal = scal*1d-2
165 :
166 0 : dudt_actual_ad = 0d0
167 0 : dudt_actual_ad%val = s% dxh_u(k)/dt
168 0 : dudt_actual_ad%d1Array(i_v_00) = 1d0/dt
169 :
170 0 : resid_ad = scal*(dudt_expected_ad - dudt_actual_ad)
171 0 : residual = resid_ad%val
172 0 : s% equ(i_du_dt, k) = residual
173 :
174 0 : if (is_bad(residual)) then
175 0 : ierr = -1
176 0 : return
177 : !$omp critical (dudt_eqn)
178 : write(*,2) 'residual', k, residual
179 : call mesa_error(__FILE__,__LINE__,'do1_dudt_eqn')
180 : !$omp end critical (dudt_eqn)
181 : end if
182 :
183 0 : call save_eqn_residual_info(s, k, nvar, i_du_dt, resid_ad, 'do1_dudt_eqn', ierr)
184 :
185 : if (test_partials) then
186 : s% solver_test_partials_val = resid_ad% val
187 : end if
188 :
189 0 : if (test_partials) then
190 : s% solver_test_partials_var = s% i_lnR
191 : s% solver_test_partials_dval_dx = resid_ad% d1Array(i_lnR_00)
192 : !write(*,*) 'do1_dudt_eqn', s% solver_test_partials_var
193 : end if
194 :
195 : contains
196 :
197 0 : subroutine setup_momentum_flux
198 0 : if (k == 1) then
199 0 : flux_out_ad = P_surf_ad*area_00
200 : else
201 0 : flux_out_ad = s% P_face_ad(k)*area_00
202 : end if
203 0 : if (k < nz) then
204 0 : flux_in_ad = shift_p1(s% P_face_ad(k+1))*area_p1
205 : else
206 0 : flux_in_ad = 0d0
207 : end if
208 0 : end subroutine setup_momentum_flux
209 :
210 0 : subroutine setup_geometry_source(ierr)
211 : use star_utils, only: calc_Ptot_ad_tw
212 : integer, intent(out) :: ierr
213 : type(auto_diff_real_star_order1) :: P
214 0 : real(dp), dimension(s% species) :: d_Ptot_dxa
215 : logical, parameter :: skip_Peos = .false., skip_mlt_Pturb = .false.
216 : ierr = 0
217 : ! use same P here as the cell pressure in P_face calculation
218 0 : call calc_Ptot_ad_tw(s, k, skip_Peos, skip_mlt_Pturb, P, d_Ptot_dxa, ierr)
219 0 : if (ierr /= 0) return
220 0 : if (k == nz) then
221 : ! no flux in from left, so only have geometry source on right
222 : ! this matters for cases with R_center > 0.
223 0 : geometry_source_ad = P*area_00
224 : else
225 0 : geometry_source_ad = P*(area_00 - area_p1)
226 : end if
227 0 : end subroutine setup_geometry_source
228 :
229 0 : subroutine setup_gravity_source
230 : type(auto_diff_real_star_order1) :: G00, Gp1, gsL, gsR
231 : real(dp) :: mR, mL
232 : ! left 1/2 of dm gets gravity force at left face
233 : ! right 1/2 of dm gets gravity force at right face.
234 : ! this form is to match the gravity force equilibrium reconstruction.
235 0 : mR = s% m(k)
236 0 : if (k == nz) then
237 0 : mL = s% M_center
238 : else
239 0 : mL = s% m(k+1)
240 : end if
241 0 : call get_G(s, k, G00)
242 0 : gsR = -G00*mR*0.5d0*dm*inv_R2_00
243 0 : if (k == nz) then
244 0 : gsL = 0d0
245 : else
246 0 : call get_G(s, k+1, Gp1)
247 0 : Gp1 = shift_p1(Gp1)
248 0 : gsL = -Gp1*mL*0.5d0*dm*inv_R2_p1
249 : end if
250 0 : gravity_source_ad = gsL + gsR ! total gravitational force on cell
251 :
252 :
253 0 : end subroutine setup_gravity_source
254 :
255 0 : subroutine setup_diffusion_source
256 : type(auto_diff_real_star_order1) :: u_m1, u_00, u_p1
257 0 : real(dp) :: sig00, sigp1
258 0 : do_diffusion = s% RTI_flag .and. s% dudt_RTI_diffusion_factor > 0d0
259 0 : if (do_diffusion) then ! add diffusion source term to dudt
260 0 : u_p1 = 0d0 ! sets val and d1Array to 0
261 0 : if (k < nz) then
262 0 : sigp1 = s% dudt_RTI_diffusion_factor*s% sig_RTI(k+1)
263 0 : u_p1%val = s% u(k+1)
264 0 : u_p1%d1Array(i_v_p1) = 1d0
265 : else
266 0 : sigp1 = 0
267 : end if
268 0 : u_m1 = 0d0 ! sets val and d1Array to 0
269 0 : if (k > 1) then
270 0 : sig00 = s% dudt_RTI_diffusion_factor*s% sig_RTI(k)
271 0 : u_m1%val = s% u(k-1)
272 0 : u_m1%d1Array(i_v_m1) = 1d0
273 : else
274 0 : sig00 = 0
275 : end if
276 0 : u_00 = 0d0 ! sets val and d1Array to 0
277 0 : u_00%val = s% u(k)
278 0 : u_00%d1Array(i_v_00) = 1d0
279 0 : diffusion_source_ad = sig00*(u_m1 - u_00) - sigp1*(u_00 - u_p1)
280 : else
281 0 : diffusion_source_ad = 0d0
282 : end if
283 0 : s% dudt_RTI(k) = diffusion_source_ad%val/dm
284 0 : end subroutine setup_diffusion_source
285 :
286 : end subroutine do1_dudt_eqn
287 :
288 :
289 0 : subroutine do_uface_and_Pface(s, ierr)
290 : type (star_info), pointer :: s
291 : integer, intent(out) :: ierr
292 : integer :: k, op_err
293 : include 'formats'
294 0 : ierr = 0
295 0 : !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2)
296 : do k = 1, s% nz
297 : op_err = 0
298 : call do1_uface_and_Pface(s, k, op_err)
299 : if (op_err /= 0) ierr = op_err
300 : end do
301 : !$OMP END PARALLEL DO
302 0 : end subroutine do_uface_and_Pface
303 :
304 :
305 0 : subroutine get_G(s, k, G)
306 : type (star_info), pointer :: s
307 : integer, intent(in) :: k
308 : type(auto_diff_real_star_order1), intent(out) :: G
309 : real(dp) :: cgrav
310 0 : cgrav = s% cgrav(k)
311 0 : G = cgrav
312 0 : if (s% rotation_flag .and. s% use_gravity_rotation_correction) &
313 0 : G = G*s% fp_rot(k)
314 0 : end subroutine get_G
315 :
316 :
317 0 : subroutine do1_uface_and_Pface(s, k, ierr)
318 : use eos_def, only: i_gamma1, i_lnfree_e, i_lnPgas
319 : use star_utils, only: calc_Ptot_ad_tw, get_face_weights
320 : use hydro_rsp2, only: compute_Uq_face
321 : type (star_info), pointer :: s
322 : integer, intent(in) :: k
323 : integer, intent(out) :: ierr
324 : logical :: test_partials
325 :
326 : type(auto_diff_real_star_order1) :: &
327 : r_ad, A_ad, PL_ad, PR_ad, uL_ad, uR_ad, rhoL_ad, rhoR_ad, &
328 : gamma1L_ad, gamma1R_ad, csL_ad, csR_ad, G_ad, dPdm_grav_ad, &
329 : Sl1_ad, Sl2_ad, Sr1_ad, Sr2_ad, numerator_ad, denominator_ad, &
330 : Sl_ad, Sr_ad, Ss_ad, P_face_L_ad, P_face_R_ad, du_ad, Uq_ad
331 0 : real(dp), dimension(s% species) :: d_Ptot_dxa ! skip this
332 : logical, parameter :: skip_Peos = .false., skip_mlt_Pturb = .false.
333 0 : real(dp) :: delta_m, f
334 :
335 : include 'formats'
336 :
337 0 : ierr = 0
338 0 : test_partials = .false.
339 : !test_partials = (k == s% solver_test_partials_k)
340 :
341 0 : s% RTI_du_diffusion_kick(k) = 0d0
342 0 : s% d_uface_domega(k) = 0
343 :
344 0 : if (k == 1) then
345 0 : s% u_face_ad(k) = wrap_u_00(s,k)
346 0 : s% P_face_ad(k) = wrap_Peos_00(s,k)
347 0 : return
348 : end if
349 :
350 0 : r_ad = wrap_r_00(s,k)
351 0 : A_ad = 4d0*pi*pow2(r_ad)
352 :
353 0 : call calc_Ptot_ad_tw(s, k, skip_Peos, skip_mlt_Pturb, PL_ad, d_Ptot_dxa, ierr)
354 0 : if (ierr /= 0) return
355 0 : call calc_Ptot_ad_tw(s, k-1, skip_Peos, skip_mlt_Pturb, PR_ad, d_Ptot_dxa, ierr)
356 0 : if (ierr /= 0) return
357 0 : PR_ad = shift_m1(PR_ad)
358 :
359 0 : uL_ad = wrap_u_00(s,k)
360 0 : uR_ad = wrap_u_m1(s,k)
361 :
362 0 : rhoL_ad = wrap_d_00(s,k)
363 0 : rhoR_ad = wrap_d_m1(s,k)
364 :
365 0 : gamma1L_ad = wrap_gamma1_00(s,k)
366 0 : gamma1R_ad = wrap_gamma1_m1(s,k)
367 :
368 0 : csL_ad = sqrt(gamma1L_ad*PL_ad/rhoL_ad)
369 0 : csR_ad = sqrt(gamma1R_ad*PR_ad/rhoR_ad)
370 :
371 : ! change PR and PL for gravity
372 0 : call get_G(s, k, G_ad)
373 :
374 0 : dPdm_grav_ad = -G_ad*s% m_grav(k)/(pow2(r_ad)*A_ad) ! cm^-1 s^-2
375 :
376 0 : delta_m = 0.5d0*s% dm(k) ! positive delta_m from left center to edge
377 0 : PL_ad = PL_ad + delta_m*dPdm_grav_ad
378 :
379 0 : delta_m = -0.5d0*s% dm(k-1) ! negative delta_m from right center to edge
380 0 : PR_ad = PR_ad + delta_m*dPdm_grav_ad
381 :
382 : ! acoustic wavespeeds (eqn 2.38)
383 0 : Sl1_ad = uL_ad - csL_ad
384 0 : Sl2_ad = uR_ad - csR_ad
385 :
386 : ! take Sl = min(Sl1, Sl2)
387 0 : if (Sl1_ad%val < Sl2_ad%val) then
388 0 : Sl_ad = Sl1_ad
389 : else
390 0 : Sl_ad = Sl2_ad
391 : end if
392 :
393 0 : Sr1_ad = uR_ad + csR_ad
394 0 : Sr2_ad = uL_ad + csL_ad
395 :
396 : ! take Sr = max(Sr1, Sr2)
397 0 : if (Sr1_ad%val > Sr2_ad%val) then
398 0 : Sr_ad = Sr1_ad
399 : else
400 0 : Sr_ad = Sr2_ad
401 : end if
402 :
403 : ! contact velocity (eqn 2.20)
404 0 : numerator_ad = uR_ad*rhoR_ad*(Sr_ad - uR_ad) + uL_ad*rhoL_ad*(uL_ad - Sl_ad) + (PL_ad - PR_ad)
405 0 : denominator_ad = rhoR_ad*(Sr_ad - uR_ad) + rhoL_ad*(uL_ad - Sl_ad)
406 :
407 0 : if (denominator_ad%val == 0d0 .or. is_bad(denominator_ad%val)) then
408 0 : ierr = -1
409 0 : if (s% report_ierr) then
410 0 : write(*,2) 'u_face denominator bad', k, denominator_ad%val
411 : end if
412 0 : return
413 : end if
414 :
415 0 : Ss_ad = numerator_ad/denominator_ad
416 :
417 0 : s% u_face_ad(k) = Ss_ad
418 0 : s% d_uface_domega(k) = s% u_face_ad(k)%d1Array(i_L_00)
419 :
420 : ! contact pressure (eqn 2.19)
421 0 : P_face_L_ad = rhoL_ad*(uL_ad-Sl_ad)*(uL_ad-Ss_ad) + PL_ad
422 0 : P_face_R_ad = rhoR_ad*(uR_ad-Sr_ad)*(uR_ad-Ss_ad) + PR_ad
423 :
424 0 : s% P_face_ad(k) = 0.5d0*(P_face_L_ad + P_face_R_ad) ! these are ideally equal
425 :
426 0 : if (k < s% nz .and. s% RTI_flag) then
427 : if (s% eta_RTI(k) > 0d0 .and. &
428 0 : s% dlnddt_RTI_diffusion_factor > 0d0 .and. s% dt > 0d0) then
429 0 : f = s% dlnddt_RTI_diffusion_factor*s% eta_RTI(k)/s% dm_bar(k)
430 0 : du_ad = f*A_ad*(rhoL_ad - rhoR_ad) ! bump uface in direction of lower density
431 0 : s% RTI_du_diffusion_kick(k) = du_ad%val
432 0 : s% u_face_ad(k) = s% u_face_ad(k) + du_ad
433 : end if
434 : end if
435 :
436 :
437 0 : if (s% RSP2_flag) then ! include Uq in u_face, To do: implement in sources instead ~ EbF
438 0 : Uq_ad = compute_Uq_face(s, k, ierr)
439 0 : if (ierr /= 0) return
440 0 : s% u_face_ad(k) = s% u_face_ad(k) + Uq_ad
441 : end if
442 :
443 0 : s% u_face_val(k) = s% u_face_ad(k)%val
444 :
445 0 : if (s% P_face_start(k) < 0d0) then
446 0 : s% u_face_start(k) = s% u_face_val(k)
447 0 : s% P_face_start(k) = s% P_face_ad(k)%val
448 : end if
449 :
450 : if (test_partials) then
451 : s% solver_test_partials_val = PL_ad% val
452 : s% solver_test_partials_var = s% i_w_div_wc
453 : s% solver_test_partials_dval_dx = PL_ad% d1Array(i_w_div_wc_00)
454 : write(*,*) 'do1_uface_and_Pface', s% solver_test_partials_var, PL_ad% val
455 : end if
456 :
457 0 : end subroutine do1_uface_and_Pface
458 :
459 : end module hydro_riemann
|
