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