Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010 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 eps_mdot
21 :
22 : use star_private_def
23 : use const_def, only: dp
24 : use star_utils
25 : use accurate_sum ! Provides the accurate_real type, which enables us to do
26 : !sums and differences without much loss of precision.
27 : use mass_utils ! Helper methods for accurately computiong mass differences
28 : ! and intersections between different meshes.
29 :
30 : implicit none
31 :
32 : private
33 : public :: calculate_eps_mdot
34 :
35 : contains
36 :
37 : !> We choose as a convention to fix F on the faces between cells and to be
38 : !! positive when mass is flowing downward. With this, the mass flux may be
39 : !! related to the change in the mass of the cell as
40 : !!
41 : !! Delta(dm_j) = F_{j} - F_{j+1}.
42 : !!
43 : !! The left-hand side is known: we have
44 : !!
45 : !! Delta(dm_j) = (dm_j_after_adust_mass - dm_j_before_adjust_mass)
46 : !!
47 : !! This is what we call change_in_dm.
48 : !!
49 : !! Furthermore we fix the flux through the centre (F_{nz+1}) to zero.
50 : !! Hence
51 : !!
52 : !! F_{nz} = Delta(dm_nz)
53 : !! F_{nz-1} = F_{nz} + Delta(dm_{nz-1})
54 : !!
55 : !! and so on.
56 : !!
57 : !! With a little rearranging we therefore find_mass_flux
58 : !!
59 : !! F_{j} = F_{j+1} + (dm_j_after_adust_mass - dm_j_before_adjust_mass)
60 : !!
61 : !! @param change_in_dm Change in cell masses.
62 : !! @param nz Number of cells.
63 : !! @param mass_flux Mass flux between cells.
64 0 : subroutine find_mass_flux(nz, change_in_dm, mass_flux)
65 : ! Inputs
66 : real(qp), dimension(:), intent(in) :: change_in_dm
67 : integer, intent(in) :: nz
68 :
69 : ! Intermediates
70 : integer :: j
71 :
72 : real(qp), dimension(:), intent(out) :: mass_flux
73 :
74 0 : mass_flux(nz+1) = 0 ! Fixes the inner boundary.
75 :
76 0 : do j=nz,1,-1
77 0 : mass_flux(j) = mass_flux(j+1) + change_in_dm(j)
78 : end do
79 :
80 0 : end subroutine find_mass_flux
81 :
82 :
83 0 : real(dp) function interpolate_onto_faces(vec, dm, nz, j)
84 : ! Inputs
85 : real(dp), dimension(:) :: vec
86 : real(qp), dimension(:) :: dm
87 : integer :: nz, j
88 :
89 : ! Intermediates
90 0 : real(dp) :: alpha, beta, dmbar1
91 :
92 : ! Return
93 : real(dp) :: interp
94 :
95 : !!! High-level explanation
96 :
97 : ! Returns the value of vec interpolated onto face j.
98 : ! This works for j running from 1 through length(vec) + 1 == nz + 1,
99 : ! so that we include the 'bounding' faces.
100 :
101 : !!! Mechanics explanation
102 :
103 : ! The interpolation is done so that the finite difference
104 : !
105 : ! (interp_vec(j+1)-interp_vec(j) / dm(j))
106 : !
107 : ! Equals the average of the finite differences
108 : !
109 : ! (vec(j+1)-vec(j))/(dm-bar(j+1))
110 : !
111 : ! and
112 : !
113 : ! (vec(j)-vec(j-1))/(dm-bar(j)),
114 : !
115 : ! where
116 : !
117 : ! dm-bar(j) = (1/2)(dm(j-1) + dm(j)).
118 : !
119 : ! This is done because these finite differences are what
120 : ! MESA is using elsewhere, so in order to ensure consistency
121 : ! we want our interpolated vector to have derivatives consistent
122 : ! with these.
123 : !
124 : ! When j == 1 we can't do this because we don't know vec(j-1), so
125 : ! we take vec(j-1) == vec(j) and return vec(1).
126 : ! When j == length(vec) + 1 we likewise need to assume vec(j-1) == vec(j)
127 : ! so we return vec(nz).
128 :
129 0 : interp = 0
130 0 : if (j == 1) interp = vec(1)
131 0 : if (j == nz + 1) interp = vec(nz)
132 0 : if (j > 1 .and. j < nz + 1) then
133 0 : dmbar1 = (dm(j-1) + dm(j)) / 2.0d0
134 0 : beta = dm(j-1) / (2.0d0 * dmbar1)
135 0 : alpha = 1 - beta
136 0 : interp = alpha * vec(j-1) + beta * vec(j)
137 : end if
138 0 : interpolate_onto_faces = interp
139 :
140 0 : end function interpolate_onto_faces
141 :
142 0 : subroutine set_leak_frac(nz, L, dm, dt, thermal_energy, grad_r_sub_grad_a, mass_flux, leak_frac, eps_mdot_leak_frac_factor)
143 : ! Inputs
144 : real(qp), dimension(:) :: mass_flux, dm
145 : real(dp), dimension(:) :: L, grad_r_sub_grad_a, thermal_energy, leak_frac
146 : real(dp) :: eps_mdot_leak_frac_factor
147 : real(dp) :: dt
148 : integer :: nz
149 :
150 : ! Intermediates
151 : integer :: k, km1
152 0 : real(dp) :: mass_flux_bar
153 :
154 : !!! High-level explanation
155 :
156 : ! This subroutine calculates the maximum fractional change in entropy that thermal
157 : ! leakage can produce. When this is zero all changes in state must be adiabatic.
158 : ! When this is one there is enough time for arbitrarily large changes in entropy.
159 :
160 : ! The derivation is as follows.
161 : ! The perturbation to the luminosity L of a fluid element due to a perturbation
162 : ! in its specific entropy S follows dL/L ~ d(grad S)/grad S. In radiative zones this
163 : ! is ~dS/S. In convective zones it's ~ dS / (S * grad_r_sub_grad_ad).
164 : ! Because grad_r_sub_grad_ad is of order unity in radiative zones we can use this correction everywhere.
165 : ! The takes to leak out is therefore dm * TdS/dL ~ dm * grad_r_sub_grad_ad * TS/L, where dm is the mass of the element.
166 : ! This is the thermal time-scale of the fluid element. By comparison, in time dm*dt/mass_flux the
167 : ! element will no longer be in the same cell, so while in this cell a fraction L*dt/(TS mass_flux)
168 : ! of this perturbation can leak out. The fraction which leaks out cannot exceed one, and so
169 : ! we obtain the expression below. Note that we approximate TS by the specific thermal energy Cp*T.
170 :
171 0 : do k=1,nz
172 0 : if (k == 1) then
173 0 : mass_flux_bar = mass_flux(1)
174 : else
175 0 : km1 = k-1
176 0 : mass_flux_bar = (mass_flux(km1) + mass_flux(k)) / 2.0d0
177 : end if
178 0 : if (abs(grad_r_sub_grad_a(k)) == 0d0 .or. mass_flux_bar == 0d0) then
179 0 : leak_frac(k) = 1.0d0
180 : else
181 : leak_frac(k) = eps_mdot_leak_frac_factor *&
182 0 : abs(L(k) * dt / (grad_r_sub_grad_a(k) * thermal_energy(k) * mass_flux_bar))
183 0 : leak_frac(k) = min(1.0d0, leak_frac(k))
184 : end if
185 : end do
186 :
187 0 : end subroutine set_leak_frac
188 :
189 0 : subroutine leak_control(nz, mass_flux, dm, mesh_intersects, ranges,&
190 0 : total_mass_through_cell, eps_mdot_per_total_mass,&
191 0 : accumulated, mdot_adiabatic_surface, leak_frac)
192 : ! Inputs
193 : integer :: nz
194 : integer, dimension(:,:) :: ranges
195 : real(dp) :: mdot_adiabatic_surface
196 : real(qp), dimension(:) :: mass_flux, dm, mesh_intersects, total_mass_through_cell
197 : real(dp), dimension(:) :: eps_mdot_per_total_mass, accumulated, leak_frac
198 :
199 : ! Intermediates
200 : integer :: j
201 : integer :: i_start, i_end
202 0 : integer, dimension(:), allocatable :: i_min, i_max, j_min, j_max
203 : real(qp) :: delta_m
204 0 : real(dp) :: sgn
205 0 : type(non_rect_array), dimension(:), allocatable :: pf
206 :
207 : !!! High-level explanation
208 :
209 : ! When the thermal time-scale is long, material changes adiabatically as it descends.
210 : ! When the thermal time-scale is short, material bleeds entropy along the way.
211 : ! The adjust_mass routine assumes that we are in the latter limit and so produces a state
212 : ! which assumes entropy has been lost along the way.
213 : ! If this is correct, we should count up the entropy lost in each cell and use that to set eps_mdot.
214 : ! If this is not correct, we should instead follow parcels of material as they descend, track
215 : ! the entropy they were (incorrectly) assumed to have lost, and restore it via eps_mdot.
216 : ! We connect these limits using the leak fraction, as defined in set_leak_frac.
217 :
218 : ! The calculation of what energy to restore and what to leak is done in the leak subroutine.
219 : ! This (leak_control) subroutine precomputes some useful quantities for that calculation and
220 : ! then divides the star into contiguous blocks within which the mass flows monotonically
221 : ! up (surface-ward) or down (center-ward). It then calls leak on each such block.
222 :
223 : ! Initialize
224 0 : accumulated = 0
225 0 : mdot_adiabatic_surface = 0
226 0 : delta_m = mass_flux(1)
227 0 : if (mass_flux(1) == 0) then
228 : ! Should never happen, because we bail out earlier if mdot == 0.
229 : ! Nevertheless the logic that follows should work without issue.
230 : ! We just need to set delta_m to be something non-zero for the purposes
231 : ! of this subroutine because -delta_m is used to normalize the pass fraction.
232 : ! That is ALL it is used to do in this routine however, and the normalization
233 : ! is arbitrary (we just choose to use delta_m for convenience) so we can set it
234 : ! to anything we like.
235 0 : delta_m = -1.
236 : end if
237 :
238 : ! Calculate pass fraction:
239 : ! pf holds the fraction of mass in cell j which passed through cell i for
240 : ! all (i,j) for which this is neither 0 nor 1. The method compute_pass_fraction
241 : ! takes this as an input and does the relevant tests to determine which of
242 : ! (0, 1, pf(j)%arr(i)) is the appropriate fraction.
243 : ! Note that the pass fraction for cell j == 0 is either zero (if delta_m > 0) or
244 : ! else is normalized against the total mass which leaves the model (-delta_m == -mass_flux(1)).
245 0 : allocate(i_min(0:nz), i_max(0:nz), pf(0:nz))
246 0 : call prepare_pass_fraction(nz, delta_m, dm, mesh_intersects, ranges, i_min, i_max, pf)
247 :
248 :
249 : ! We determine for each cell i the set of cells {j} whose ending material
250 : ! passed through i during adjust_mass. This is a contiguous set because if
251 : ! material in cells k1 and k2 (k1 < k2) passed through cell i then because
252 : ! the material between cells k1 and k2 is between them it must also have passed
253 : ! through cell i. Hence we just need to find the first and last cell which
254 : ! contains material which passed through cell i. j_min and j_max. These are
255 : ! computed in the mass_utils subroutine find_j_ranges.
256 0 : allocate(j_min(nz), j_max(nz))
257 0 : call find_j_ranges(nz, ranges, mass_flux, j_min, j_max)
258 :
259 : ! Starting state
260 0 : i_start = 1 ! Where we'll begin leaking
261 0 : i_end = 1 ! Where we'll end leaking
262 0 : if (mass_flux(1) == 0) then
263 : sgn = 0
264 : else
265 0 : sgn = mass_flux(1) / abs(mass_flux(1))
266 : end if
267 :
268 : ! The following loop cuts the star up into segments in which
269 : ! mass is either all flowing or all flowing down. It then calls
270 : ! the leak subroutine on each of these in turn.
271 0 : do while (i_start <= nz .and. i_end <= nz)
272 0 : if (sgn > 0) then
273 0 : if (mass_flux(i_end + 1) > 0) then
274 0 : i_end = i_end + 1
275 0 : else if (mass_flux(i_end + 1) < 0) then
276 : call leak(nz, i_start, i_end, i_min, i_max, j_min, j_max, pf,&
277 : dm, delta_m, accumulated, eps_mdot_per_total_mass, leak_frac,&
278 0 : total_mass_through_cell, mdot_adiabatic_surface)
279 0 : i_start = i_end
280 0 : sgn = -1
281 : else
282 : call leak(nz, i_start, i_end, i_min, i_max, j_min, j_max, pf,&
283 : dm, delta_m, accumulated, eps_mdot_per_total_mass, leak_frac,&
284 0 : total_mass_through_cell, mdot_adiabatic_surface)
285 0 : i_start = i_end
286 0 : sgn = 0
287 : end if
288 0 : else if (sgn < 0) then
289 0 : if (mass_flux(i_start + 1) < 0) then
290 0 : i_start = i_start + 1
291 0 : else if (mass_flux(i_start + 1) > 0) then
292 : call leak(nz, i_start, i_end, i_min, i_max, j_min, j_max, pf,&
293 : dm, delta_m, accumulated, eps_mdot_per_total_mass, leak_frac,&
294 0 : total_mass_through_cell, mdot_adiabatic_surface)
295 0 : i_end = i_start
296 0 : sgn = 1
297 : else
298 : call leak(nz, i_start, i_end, i_min, i_max, j_min, j_max, pf,&
299 : dm, delta_m, accumulated, eps_mdot_per_total_mass, leak_frac,&
300 0 : total_mass_through_cell, mdot_adiabatic_surface)
301 0 : i_end = i_start
302 0 : sgn = 0
303 : end if
304 : else
305 0 : if (mass_flux(i_start + 1) > 0) then
306 0 : i_start = i_start + 1
307 0 : i_end = i_end + 1
308 0 : sgn = 1
309 0 : else if (mass_flux(i_start + 1) < 0) then
310 0 : i_start = i_start + 1
311 0 : i_end = i_end + 1
312 0 : sgn = -1
313 : else
314 : call leak(nz, i_start, i_end, i_min, i_max, j_min, j_max, pf,&
315 : dm, delta_m, accumulated, eps_mdot_per_total_mass, leak_frac,&
316 0 : total_mass_through_cell, mdot_adiabatic_surface)
317 0 : i_start = i_start + 1
318 0 : i_end = i_end + 1
319 : end if
320 : end if
321 : end do
322 :
323 0 : do j=1,nz
324 0 : accumulated(j) = accumulated(j) / dm(j)
325 : end do
326 :
327 : ! This handles the contribution of material which starts and ends in the same cell,
328 : ! or which starts outside the star and ends in the top cell.
329 0 : do j=1,nz
330 0 : if (leak_frac(j) == 1) then
331 : ! Means all energy leaks into this cell, so we don't do the normal leak calculation.
332 0 : accumulated(j) = accumulated(j) + eps_mdot_per_total_mass(j) * total_mass_through_cell(j)
333 : else
334 : ! Leakage happens so we just need this piece.
335 0 : accumulated(j) = accumulated(j) + eps_mdot_per_total_mass(j) * dm(j) * compute_pass_fraction(j, j, i_min, i_max, pf)
336 : end if
337 : end do
338 :
339 0 : end subroutine leak_control
340 :
341 :
342 0 : subroutine leak(nz, i_start, i_end, i_min, i_max, j_min, j_max, pf,&
343 0 : dm, delta_m, accumulated, eps_mdot_per_total_mass, leak_frac,&
344 : total_mass_through_cell, mdot_adiabatic_surface)
345 :
346 : ! Inputs
347 : integer :: nz, i_start, i_end
348 : integer, dimension(:) :: i_min, i_max, j_min, j_max
349 : type(non_rect_array), dimension(:) :: pf
350 : real(qp) :: delta_m
351 : real(dp), dimension(:) :: leak_frac, accumulated, eps_mdot_per_total_mass
352 : real(qp), dimension(:) :: dm, total_mass_through_cell
353 : real(dp) :: mdot_adiabatic_surface
354 :
355 :
356 : ! Intermediates
357 : integer :: i, j, direction
358 0 : real(qp) :: pass_frac, next, pass_mass
359 0 : real(dp), dimension(:), allocatable :: excess
360 :
361 : !!! High-level explanation
362 :
363 : ! The leak fraction determines the fraction of any attempted entropy
364 : ! change which leaks out (versus being held adiabatically).
365 : ! excess is the amount of heat the material has which it would deposit
366 : ! if given the time.
367 : ! Hence if we follow the material that ultimately ends in cell j as it
368 : ! passes through cell i we find
369 : !
370 : ! excess(j) -> excess(j) + (heat picked up while passing through cell i)
371 : ! (deposited heat) = leak_frac(i) * excess(j)
372 : ! excess(j) -> (1 - leak_frac) * excess(j)
373 : !
374 : ! That is, we first calculate the entropy change associated with moving the material
375 : ! through the next cell (heat picked up). We then add that to the current excess.
376 : ! A fraction leak_frac of the cumulative entropy excess is then deposited as heat and
377 : ! decremented from the excess.
378 : ! When the material reaches whatever cell it ends in (i == j) the excess is deposited
379 : ! in that cell. If material exits the star the excess it leaves with is accounted for
380 : ! in mdot_adiabatic_surface.
381 :
382 :
383 :
384 : ! There are only two ways for i_start to equal i_end:
385 : ! 1. All mass which begins in this cell ends in this cell so there's nothing to do.
386 : ! 2. This is the top cell and all mass which ends in this cell begins in either this cell or outside the star.
387 : ! Because no eps_mdot is accumulated outside the star this is equivalent to it all
388 : ! beginning in the top cell, so there's again nothing to do.
389 : ! In either case the loop at the end of leak_control will ensure that
390 : ! accumulated(i_start) = eps_mdot_per_total_mass(i) * dm(i),
391 : ! which is just the previously calculated eps_mdot(i).
392 0 : if (i_start == i_end) return
393 :
394 : ! Determine flow direction
395 0 : if (i_start < i_end) then
396 : direction = 1
397 : else
398 0 : direction = -1
399 : end if
400 :
401 0 : allocate(excess(0:nz))
402 0 : excess = 0
403 0 : i = i_start
404 0 : do while (i /= i_end + direction)
405 0 : if (leak_frac(i) == 1) then
406 : ! If this is the first cell with full leakage we dump all excess here.
407 : ! If not it means we've done this before (so there is no excess)
408 : ! or we're at the very start (so there is no excess).
409 0 : if (i /= i_start) then
410 0 : if (leak_frac(i - direction) < 1) then
411 0 : do j=j_min(i), j_max(i)
412 0 : accumulated(i) = accumulated(i) + excess(j)
413 0 : excess(j) = 0
414 : end do
415 : end if
416 : end if
417 : else
418 0 : do j=j_min(i), j_max(i)
419 0 : if ((i == i_end .and. j > 0) .or. j == i) then
420 : ! This material ends here, so it drops its heat here.
421 : ! The first condition means that this is the last cell (i)
422 : ! to be considered in this flow direction and that cell (j)
423 : ! does not represent material that exits the star.
424 : ! The second condition means the cell under consideration (i)
425 : ! is the same as that of the material under consideration (j)
426 : ! and so that material must end up here.
427 :
428 : ! Note that we do not include the eps_mdot contribution from
429 : ! this cell here because that could cause double counting in
430 : ! cells across which mass_flux changes sign. To ensure single
431 : ! counting that contribution is accounted for in the loop
432 : ! at the end of leak_frac.
433 0 : accumulated(i) = accumulated(i) + excess(j)
434 0 : excess(j) = 0
435 0 : else if (i == i_end .and. i == 1 .and. j == 0) then
436 : ! Material with j == 0 exits the star. Note that this implies direction == -1.
437 : ! For i > 1 this material can be handled by the 'just passing through' else
438 : ! clause, so we only need to think about the i == 1 case.
439 :
440 : ! Because this material isn't in the star at the end, we have to account
441 : ! for the heat it picks up on its way out:
442 0 : pass_frac = compute_pass_fraction(1, 0, i_min, i_max, pf)
443 0 : excess(j) = excess(j) - delta_m * pass_frac * eps_mdot_per_total_mass(1)
444 :
445 : ! Because our accounting assumes that material flows through the
446 : ! surface of the star with energy equal to the surface specific energy
447 : ! we have to account for any excess this material carries on its way out.
448 0 : mdot_adiabatic_surface = mdot_adiabatic_surface + (1 - leak_frac(i_end)) * excess(j)
449 0 : accumulated(i) = accumulated(i) + leak_frac(i_end) * excess(j)
450 :
451 0 : excess(j) = 0
452 : else
453 : ! The material in cell j is just passing through cell i.
454 0 : pass_frac = compute_pass_fraction(i, j, i_min, i_max, pf)
455 0 : if (j > 0) then
456 0 : pass_mass = pass_frac * dm(j)
457 : else
458 0 : pass_mass = -pass_frac * delta_m
459 : end if
460 0 : next = excess(j) + dm(i) * eps_mdot_per_total_mass(i) * pass_mass
461 :
462 0 : accumulated(i) = accumulated(i) + leak_frac(i) * next
463 0 : excess(j) = (1 - leak_frac(i)) * next
464 : end if
465 : end do
466 : end if
467 0 : i = i + direction
468 : end do
469 :
470 0 : end subroutine leak
471 :
472 11 : subroutine calculate_eps_mdot(s, dt, ierr)
473 : use adjust_mass, only: compute_prev_mesh_dm
474 :
475 : ! Inputs
476 : type (star_info), pointer :: s
477 : real(dp) :: dt
478 : integer :: ierr
479 :
480 : ! Intermediates
481 : logical, parameter :: dbg = .false.
482 : integer :: nz, j
483 11 : real(dp) :: delta_m, change_sum, leak_sum, err, abs_err, mdot_adiabatic_surface, gradT_mid
484 : real(dp), dimension(:), allocatable :: &
485 11 : p_bar, rho_bar, te_bar, te, &
486 11 : leak_frac, thermal_energy, density_weighted_flux, eps_mdot_per_total_mass,&
487 11 : accumulated, grad_r_sub_grad_a
488 11 : real(qp), dimension(:), allocatable :: change_in_dm, mass_flux, dm, prev_mesh_dm,&
489 11 : total_mass_through_cell
490 : type(accurate_real) :: sum
491 11 : integer, dimension(:,:), allocatable :: ranges
492 11 : real(qp), dimension(:), allocatable :: mesh_intersects
493 :
494 11 : if (s% mstar_dot == 0d0 .or. dt <= 0d0) then
495 13098 : s% eps_mdot(1:s%nz) = 0d0
496 11 : s% mdot_adiabatic_surface = 0d0
497 11 : s% mdot_acoustic_surface = 0d0
498 11 : s% total_internal_energy_after_adjust_mass = 0d0
499 11 : s% total_gravitational_energy_after_adjust_mass = 0d0
500 11 : s% total_radial_kinetic_energy_after_adjust_mass = 0d0
501 11 : s% total_turbulent_energy_after_adjust_mass = 0d0
502 11 : s% total_rotational_kinetic_energy_after_adjust_mass = 0d0
503 11 : s% total_energy_after_adjust_mass = 0d0
504 11 : return
505 : end if
506 :
507 0 : s% need_to_setvars = .true.
508 :
509 : ! Stellar properties
510 0 : nz = s%nz
511 0 : delta_m = s%mstar_dot * dt
512 :
513 0 : allocate(prev_mesh_dm(nz), dm(nz), change_in_dm(nz))
514 0 : call compute_prev_mesh_dm(s, prev_mesh_dm, dm, change_in_dm)
515 :
516 0 : allocate(mass_flux(nz+1))
517 0 : call find_mass_flux(nz, change_in_dm, mass_flux)
518 :
519 : ! Tabulate cell intersection widths between the new mesh and the old
520 0 : allocate(mesh_intersects(2*nz))
521 0 : allocate(ranges(2*nz,2))
522 0 : call make_compressed_intersect(dm, prev_mesh_dm, nz, mesh_intersects, ranges)
523 :
524 :
525 : ! Weighted dm/dt by rho, used for calculating PdV work.
526 0 : allocate(density_weighted_flux(nz+1))
527 0 : density_weighted_flux(nz+1) = 0
528 0 : do j=nz,1,-1
529 0 : density_weighted_flux(j) = density_weighted_flux(j+1) + change_in_dm(j) / s%rho(j)
530 : end do
531 :
532 : ! We attribute eps_mdot evenly to all of the mass which is at any point within a cell.
533 : ! This is a zeroth-order accurate scheme, but it is good enough for our purposes.
534 : ! The next best approximation would be to track the amount of "time" material spends
535 : ! in each cell and assign contributions proportional to that, but unless there is evidence
536 : ! suggesting this is necessary it seems overly complicated.
537 0 : allocate(total_mass_through_cell(nz))
538 0 : !$OMP PARALLEL DO
539 : do j=1,nz
540 : ! This equals Sum_j pass_fraction(i,j) * dm(j)
541 : total_mass_through_cell(j) = prev_mesh_dm(j) &
542 : + max(mass_flux(j),0.0_qp) - min(mass_flux(j+1),0.0_qp)
543 : end do
544 : !$OMP END PARALLEL DO
545 :
546 : ! Puts total_energy in specific terms and interpolates that as well as p and rho onto faces.
547 : ! For convenience we include the bottom face (face nz+1) and the top face (face 1), for which the
548 : ! interpolation just returns the value in the adjacent cell.
549 0 : allocate(te(nz))
550 0 : allocate(p_bar(nz+1))
551 0 : allocate(rho_bar(nz+1))
552 0 : allocate(te_bar(nz+1))
553 :
554 0 : !$OMP PARALLEL DO
555 : do j=1,nz
556 : te(j) = s% total_energy_profile_before_adjust_mass(j) / prev_mesh_dm(j)
557 : end do
558 : !$OMP END PARALLEL DO
559 0 : call eval_total_energy_profile(s, s% total_energy_profile_after_adjust_mass)
560 :
561 :
562 0 : !$OMP PARALLEL DO
563 : do j=1,nz+1
564 : ! We use the previous mesh for interpolation because that's the one for which our derivatives were calculated.
565 : p_bar(j) = interpolate_onto_faces(s%Peos, prev_mesh_dm, nz, j)
566 : rho_bar(j) = interpolate_onto_faces(s%rho, prev_mesh_dm, nz, j)
567 : te_bar(j) = interpolate_onto_faces(te, prev_mesh_dm, nz, j)
568 : end do
569 : !$OMP END PARALLEL DO
570 :
571 0 : if (s% use_other_accreting_state) then
572 0 : call s%other_accreting_state(s%id, te_bar(1), p_bar(1), rho_bar(1), ierr)
573 0 : s% surface_cell_specific_total_energy_old = te_bar(1)
574 : end if
575 :
576 : ! Calculate grad_r_sub_grad_a
577 0 : allocate(grad_r_sub_grad_a(nz))
578 0 : !$OMP PARALLEL DO PRIVATE(gradT_mid)
579 : do j=1,nz-1
580 : gradT_mid = 0.5d0*(s% gradT(j) + s% gradT(j+1))
581 : grad_r_sub_grad_a(j) = s% grada(j) - gradT_mid
582 : end do
583 : !$OMP END PARALLEL DO
584 :
585 0 : grad_r_sub_grad_a(nz) = grad_r_sub_grad_a(nz-1)
586 :
587 : ! Get thermal energy, which is not quite the same as the specific internal energy
588 : ! because there may be an offset on that (e.g. due to a high Fermi level).
589 : ! What we want, more specifically, is a quantity which, when perturbed, perturbs
590 : ! the luminosity to a comparable degree. S*T is a good approximation of this.
591 0 : allocate(thermal_energy(nz))
592 0 : !$OMP PARALLEL DO
593 : do j=1,nz
594 : thermal_energy(j) = s%T(j) * s%cp(j)
595 : end do
596 : !$OMP END PARALLEL DO
597 :
598 : ! Our PV and TE methods can only handle a situation in which all cells increase in mass.
599 : ! We handle the case in which all cells decrease in mass by swapping dm and prev_mesh_dm.
600 : ! In effect this let's us treat mass loss as mass accretion in reverse. None of the thermodynamic
601 : ! routines care about this ordering so we can just negate the output. The leakage routine
602 : ! does care, because it undergoes an irreversible process, so to make sure that is handled correctly
603 : ! we explicitly pass the direction of the change to that routine.
604 :
605 0 : allocate(eps_mdot_per_total_mass(nz))
606 0 : do j=1,nz
607 0 : eps_mdot_per_total_mass(j) = 0
608 0 : s%eps_mdot(j) = 0
609 : end do
610 :
611 0 : change_sum = 0
612 :
613 : ! Calculate change in heat.
614 0 : s % mdot_acoustic_surface = delta_m * p_bar(1) / rho_bar(1) / dt
615 0 : do j=1,nz
616 : ! We calculate eps_mdot using accurate reals to reduce roundoff errors.
617 0 : sum % sum = 0.0d0
618 0 : sum % compensator = 0.0d0
619 :
620 : ! More numerically stable version of
621 : ! te_bar(j) * mass_flux(j) - te_bar(j+1) * mass_flux(j+1) - change_in_dm(j) * te(j)
622 0 : sum = sum + real(change_in_dm(j) * te_bar(j),qp)
623 0 : sum = sum + real(mass_flux(j+1) * (te_bar(j) - te_bar(j+1)),qp)
624 0 : sum = sum - real(change_in_dm(j) * te(j),qp)
625 :
626 : ! (PV F_m - PV_0)_{k} - (PV F_m - PV_0)_{k+1}
627 0 : sum = sum + real(p_bar(j) * (mass_flux(j) / rho_bar(j) - density_weighted_flux(j)),qp)
628 0 : sum = sum - real(p_bar(j+1) * (mass_flux(j+1) / rho_bar(j+1) - density_weighted_flux(j+1)),qp)
629 :
630 : ! Change in specific total energy in the cell with the final cell mass. The
631 : ! In the absence of composition changes this is just dm_new * (change in gravitational potential),
632 : ! but when composition changes this term also captures that.
633 0 : sum = sum - real(s% total_energy_profile_after_adjust_mass(j) - s%dm(j) * te(j), qp)
634 :
635 0 : change_sum = change_sum + sum%value() / (dt)
636 :
637 : ! Multiplicative factors
638 0 : eps_mdot_per_total_mass(j) = sum % value() / (s%dm(j) * dt) / total_mass_through_cell(j)
639 :
640 : end do
641 :
642 0 : err = 0d0
643 0 : do j=1,nz
644 : err = err + eps_mdot_per_total_mass(j) * s%dm(j) * dt * total_mass_through_cell(j)
645 : end do
646 0 : err = err - s%mdot_acoustic_surface
647 0 : err = err - te_bar(1) * delta_m
648 :
649 : ! Calculate leak fraction
650 0 : allocate(leak_frac(nz))
651 : call set_leak_frac(nz, s%L, dm, dt, thermal_energy, grad_r_sub_grad_a,&
652 0 : mass_flux, leak_frac, s%eps_mdot_leak_frac_factor)
653 :
654 : ! Now we calculate what ought to have happened accounting for energy retention/leakage.
655 0 : allocate(accumulated(nz))
656 : call leak_control(nz, mass_flux, dm, mesh_intersects, ranges,&
657 : total_mass_through_cell, eps_mdot_per_total_mass,&
658 0 : accumulated, mdot_adiabatic_surface, leak_frac)
659 0 : do j=1,nz
660 0 : s%eps_mdot(j) = accumulated(j)
661 : end do
662 0 : s% mdot_adiabatic_surface = -mdot_adiabatic_surface
663 :
664 : ! Apply user-specified scaling factor.
665 : ! This non-physical, but can be useful for stellar engineering.
666 0 : s% eps_mdot(1:nz) = s% eps_mdot_factor * s% eps_mdot(1:nz)
667 :
668 : ! Evaluate total energies for later use outside of eps_mdot.
669 : call eval_deltaM_total_energy_integrals( &
670 : s, 1, nz, s% mstar, .true., &
671 : s% total_energy_profile_after_adjust_mass, &
672 : s% total_internal_energy_after_adjust_mass, &
673 : s% total_gravitational_energy_after_adjust_mass, &
674 : s% total_radial_kinetic_energy_after_adjust_mass, &
675 : s% total_rotational_kinetic_energy_after_adjust_mass, &
676 : s% total_turbulent_energy_after_adjust_mass, &
677 0 : s% total_energy_after_adjust_mass)
678 :
679 : if (dbg) then
680 : write(*,*) 'Mdot',s%mstar_dot
681 : leak_sum = mdot_adiabatic_surface
682 : do j=1,nz
683 : leak_sum = leak_sum + s%eps_mdot(j) * s%dm(j)
684 : end do
685 : write(*,*) 'Leak Difference:', (leak_sum - change_sum) / change_sum, dt * change_sum, dt * leak_sum, &
686 : dt * (leak_sum - change_sum)
687 :
688 : err = 0
689 : abs_err = 0
690 : do j=1,nz
691 : err = err + s% total_energy_profile_after_adjust_mass(j) - te(j) * prev_mesh_dm(j)
692 : err = err + s%eps_mdot(j) * dm(j) * dt
693 : ! abs_err = abs_err + abs(s% total_energy_profile_after_adjust_mass(j) - te(j) * prev_mesh_dm(j))
694 : abs_err = abs_err + abs(s% total_energy_profile_after_adjust_mass(j))
695 : ! write(*,"(5ES16.8)") s% total_energy_profile_after_adjust_mass(j), te(j), prev_mesh_dm(j), &
696 : ! s% total_energy_profile_after_adjust_mass(j) - te(j) * prev_mesh_dm(j), s%eps_mdot(j) * dm(j) * dt
697 : end do
698 :
699 : err = err - te(1) * delta_m - (s%mdot_acoustic_surface + s%mdot_adiabatic_surface) * dt
700 : write(*,*) 'Err:', err, err / abs_err, err / delta_m, delta_m / s%m(1), abs_err, s%mdot_acoustic_surface*dt
701 : end if
702 :
703 :
704 11 : end subroutine calculate_eps_mdot
705 :
706 : end module eps_mdot
|
