Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2013 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 mass_utils
21 :
22 : use const_def, only: dp, qp
23 : use accurate_sum ! Provides the accurate_real type, which enables us to do
24 : !sums and differences without much loss of precision.
25 :
26 : implicit none
27 :
28 : private
29 : public :: accurate_mass_difference
30 : public :: reconstruct_m
31 : public :: reconstruct_xm
32 : public :: make_compressed_intersect
33 : public :: non_rect_array
34 : public :: prepare_pass_fraction
35 : public :: compute_pass_fraction
36 : public :: find_j_ranges
37 : public :: integrate_conserved
38 :
39 : ! A derived array type for storing non-rectangular 2D arrays
40 : type non_rect_array
41 : real(qp), dimension(:), allocatable :: arr
42 : end type non_rect_array
43 :
44 : contains
45 :
46 : ! Reconstructs m(j) given dm.
47 : ! Not currently used, but helpful for debugging.
48 0 : real(qp) function reconstruct_m(dm, nz, j)
49 : ! Inputs
50 : real(qp), dimension(:) :: dm
51 : integer :: nz, j
52 :
53 : ! Intermediates
54 : real(qp) :: sum, compensator
55 : integer :: l
56 :
57 0 : sum = 0.0
58 0 : compensator = 0.0
59 0 : do l=nz,j,-1
60 0 : call neumaier_sum(sum, compensator, dm(l))
61 : end do
62 :
63 0 : reconstruct_m = sum + compensator
64 0 : end function reconstruct_m
65 :
66 :
67 : ! Reconstructs xm(j) given dm.
68 : ! Not currently used, but helpful for debugging.
69 0 : real(qp) function reconstruct_xm(dm, nz, j)
70 : ! Inputs
71 : real(qp), dimension(:) :: dm
72 : integer :: nz, j
73 :
74 : ! Intermediates
75 : real(qp) :: sum, compensator
76 : integer :: l
77 :
78 0 : sum = 0.0
79 0 : compensator = 0.0
80 0 : do l=1,j
81 0 : call neumaier_sum(sum, compensator, dm(l))
82 : end do
83 :
84 0 : reconstruct_xm = sum + compensator
85 0 : end function reconstruct_xm
86 :
87 :
88 : ! Returns
89 : !
90 : ! m1(j) - m2(k)
91 : !
92 : ! where
93 : !
94 : ! m_i(j) = Sum_{l >= j} dm_i(l)
95 : !
96 : ! The reason we specify dm's rather than m's is to avoid
97 : ! having to subtract large numbers and incur the resultant
98 : ! penalty in precision.
99 : !
100 : ! We require length(dm1) == length(dm2) == nz.
101 : ! We permit j and k to range from 1 ... nz+1 inclusive.
102 : ! By the above definition m_i(nz+1) == 0.
103 : !
104 : ! With the above we begin by writing
105 : !
106 : ! m1(j) - m2(k) = Sum_{l >= j} dm1(l) - Sum_{l >= k} dm2(l)
107 : !
108 : ! We then apply Neumaier's algorithm (see accurate_real docs)
109 : ! to perform this sum accurately. As a preprocessing step we present
110 : ! terms to this algorithm with alternating signs and
111 : ! descending magnitudes for as long as possible.
112 : ! Hence we actually begin with the
113 : !
114 : ! Sum_{l = nz ... max(j,k)} dm1(l) - dm2(l) ()
115 : !
116 : ! If j < k we then add Sum_{l=k-1 ... j} dm1(l).
117 : ! If j > k we then subtract Sum_{l=j-1 ... k} dm2(l).
118 : !
119 : ! Not currently used, but helpful for debugging.
120 0 : real(qp) function accurate_mass_difference(dm1, dm2, j, k, nz)
121 : ! Inputs
122 : real(qp), dimension(:) :: dm1, dm2
123 : integer :: j, k, nz
124 :
125 : ! Intermediates
126 : type(accurate_real) :: sum
127 0 : real(qp) :: summand
128 : integer :: l
129 :
130 0 : sum % sum = 0.0
131 0 : sum % compensator = 0.0
132 :
133 0 : if (max(j,k) <= nz) then
134 0 : do l=nz,max(j,k),-1
135 0 : sum = sum + dm1(l)
136 0 : sum = sum - dm2(l)
137 : end do
138 : end if
139 :
140 0 : if (j /= k) then
141 0 : do l=max(j,k) - 1,min(j,k),-1
142 0 : if (j < k) then
143 0 : summand = dm1(l)
144 : else
145 0 : summand = -dm2(l)
146 : end if
147 :
148 0 : sum = sum + summand
149 :
150 : end do
151 : end if
152 :
153 0 : accurate_mass_difference = sum % value()
154 0 : end function accurate_mass_difference
155 :
156 :
157 : ! Returns the size of the overlap between the intervals
158 : !
159 : ! [m1(j+1), m1(j)]
160 : !
161 : ! and
162 : !
163 : ! [m2(k+1), m2(k)]
164 : !
165 : ! where
166 : !
167 : ! m_i(i) = Sum_{l >= i} dm_i(l)
168 : !
169 : ! The reason we specify dm's rather than m's is to avoid
170 : ! having to subtract large numbers and incur the resultant
171 : ! penalty in precision.
172 : !
173 : ! We assume dm1 and dm2 have the same length nz.
174 : ! We allow the indices j and k to range between 0 and nz-1
175 : ! inclusive. Here m(nz) is just the mass of the centre.
176 : !
177 : ! The width of the overlap is given by
178 : !
179 : ! min(m1(j), m2(k)) - max(m1(j+1), m2(k+1))
180 : !
181 : ! if positive. Otherwise the intervals do not intersect.
182 : ! We first compute the differences between the arguments of each
183 : ! min and max to determine the final difference.
184 : !
185 : ! There are nz*nz combinations of (j,k), but only <= 2*nz of them
186 : ! are non-zero. To see this assume withouut loss of generality that
187 : ! m1(1) < m2(1). Augment m1 with m1(0) = m2(1). It follows that the
188 : ! overlap (0,1) is non-zero, as is the overlap (nz,nz). Because
189 : ! mass is monotonic, the (j,k) with non-zero overlaps form a sequence
190 : ! which is monotonic in both j and k. There are at most 2*nz points in
191 : ! such a sequence. To avoid allocating a quadratic amount of memory in nz
192 : ! we therefore store just 2*nz overlaps. The locations of these overlaps are
193 : ! given by j == ranges(:,1), k == ranges(:,2). Some of these j's or k's
194 : ! will generally be zero, corresponding to the padding described above.
195 0 : subroutine make_compressed_intersect(dm1, dm2, nz, mesh_intersects, ranges)
196 : ! Inputs
197 : real(qp), dimension(:) :: dm1, dm2
198 : real(qp), dimension(:) :: mesh_intersects
199 : integer, dimension(:,:) :: ranges
200 : integer :: nz
201 :
202 : ! Intermediates
203 : integer :: side
204 0 : real(qp), dimension(:), allocatable :: remainders1, remainders2
205 : type(accurate_real) :: diff, dBottom, dTop, dBottomTop, dTopBottom
206 : integer :: i, j, k, counter
207 :
208 0 : allocate(remainders1(nz), remainders2(nz))
209 :
210 0 : j = nz
211 0 : k = nz
212 0 : counter = 2 * nz
213 :
214 0 : mesh_intersects = 0
215 :
216 : ! j tracks our position in dm1, k tracks our position in dm2.
217 : ! We begin with both equal to nz.
218 : ! We first compute the intersection between the nz cells.
219 : ! After that our algorithm is:
220 : ! 1. Begin with whichever side has its current cell has its top face deepest down.
221 : ! Decrement that index. If it cannot be decremented
222 : ! terminate, because all remaining overlaps will be zero.
223 : ! 2. Calculate the overlap between the new cell and the old one on the other side.
224 : ! 3. Return to 1.
225 : !
226 : ! Along the way we store overlaps in the order in which we encounter them.
227 : ! As discussed in the comments above, we know there will be precisely 2*nz of these.
228 : ! The indices of these are stored as ranges(counter,1) = j, ranges(counter,2) = k.
229 0 : ranges = -1
230 :
231 : ! remainders1(j) is the amount of mass in dm1(j) above all cells of dm2.
232 : ! remainders2(k) is the amount of mass in dm2(k) above all cells of dm1.
233 : ! Only one of these may have non-zero entries, and those entries form a contiguous
234 : ! block which, if non-empty, includes the surface cell.
235 0 : remainders1 = dm1
236 0 : remainders2 = dm2
237 :
238 : ! For convenience we track
239 : ! dTop = m1(j) - m2(k)
240 : ! dBottom = m1(j+1) - m2(k+1)
241 : ! dTopBottom = m1(j) - m2(k+1)
242 : ! dBottomTop = m1(j+1) - m2(k)
243 0 : dBottom = 0.0_qp
244 0 : dTop = dm1(nz) - dm2(nz)
245 0 : dBottomTop = -dm2(nz)
246 0 : dTopBottom = dm1(nz)
247 :
248 0 : if (dm1(j) > dm2(j)) then
249 0 : side = 1
250 :
251 0 : ranges(counter, 1) = j
252 0 : ranges(counter, 2) = k
253 0 : mesh_intersects(counter) = dm2(k)
254 0 : counter = counter - 1
255 :
256 0 : remainders1(j) = remainders1(j) - dm2(k)
257 0 : remainders2(j) = remainders2(j) - dm2(k)
258 0 : diff = dm1(j) - dm2(k)
259 : else
260 0 : side = 2
261 :
262 0 : ranges(counter, 1) = j
263 0 : ranges(counter, 2) = k
264 0 : mesh_intersects(counter) = dm1(j)
265 0 : counter = counter - 1
266 :
267 0 : remainders1(j) = remainders1(j) - dm1(j)
268 0 : remainders2(j) = remainders2(j) - dm1(j)
269 0 : diff = dm2(k) - dm1(j)
270 : end if
271 :
272 0 : do while (k > 0 .and. j > 0 .and. ((side == 1 .and. k > 1) .or. (side == 2 .and. j > 1)))
273 0 : if (side == 1) then
274 0 : k = k - 1
275 0 : if (dm2(k) < diff) then
276 0 : ranges(counter, 1) = j
277 0 : ranges(counter, 2) = k
278 0 : mesh_intersects(counter) = dm2(k)
279 0 : counter = counter - 1
280 :
281 0 : remainders1(j) = remainders1(j) - dm2(k)
282 0 : remainders2(k) = remainders2(k) - dm2(k)
283 0 : diff = diff - dm2(k)
284 : else
285 0 : ranges(counter, 1) = j
286 0 : ranges(counter, 2) = k
287 0 : mesh_intersects(counter) = diff
288 0 : counter = counter - 1
289 :
290 0 : remainders1(j) = remainders1(j) - diff
291 0 : remainders2(k) = remainders2(k) - diff
292 0 : diff = dm2(k) - diff
293 0 : side = 2
294 : end if
295 :
296 0 : dBottom = dBottom - dm2(k+1)
297 0 : dTopBottom = dTopBottom - dm2(k+1)
298 0 : dTop = dTop - dm2(k)
299 0 : dBottomTop = dBottomTop - dm2(k)
300 :
301 : else
302 0 : j = j -1
303 0 : if (dm1(j) < diff) then
304 0 : ranges(counter, 1) = j
305 0 : ranges(counter, 2) = k
306 0 : mesh_intersects(counter) = dm1(j)
307 0 : counter = counter - 1
308 :
309 0 : remainders1(j) = remainders1(j) - dm1(j)
310 0 : remainders2(k) = remainders2(k) - dm1(j)
311 0 : diff = diff - dm1(j)
312 : else
313 0 : ranges(counter, 1) = j
314 0 : ranges(counter, 2) = k
315 0 : mesh_intersects(counter) = diff
316 0 : counter = counter - 1
317 :
318 0 : remainders1(j) = remainders1(j) - diff
319 0 : remainders2(k) = remainders2(k) - diff
320 0 : diff = dm1(j) - diff
321 0 : side = 1
322 : end if
323 :
324 0 : dBottom = dBottom + dm1(j+1)
325 0 : dBottomTop = dBottomTop + dm1(j+1)
326 0 : dTopBottom = dTopBottom + dm1(j)
327 0 : dTop = dTop + dm1(j)
328 :
329 : end if
330 :
331 : end do
332 :
333 0 : if (j > 1) then
334 0 : do i=j,1,-1
335 0 : ranges(counter,1) = i
336 0 : ranges(counter,2) = 0
337 0 : mesh_intersects(counter) = remainders1(i)
338 0 : counter = counter - 1
339 : end do
340 : else
341 0 : do i=k,1,-1
342 0 : ranges(counter,1) = 0
343 0 : ranges(counter,2) = i
344 0 : mesh_intersects(counter) = remainders2(i)
345 0 : counter = counter - 1
346 : end do
347 : end if
348 :
349 0 : end subroutine make_compressed_intersect
350 :
351 : ! Computes for each final cell j the first (i_min(j)) and
352 : ! last (i_max(j)) initial cells whose overlap with j is non-zero.
353 0 : subroutine find_i_ranges(nz, ranges, i_min, i_max)
354 : ! Inputs
355 : integer :: nz
356 : integer, dimension(:,:) :: ranges
357 : integer, dimension(0:nz) :: i_min, i_max
358 :
359 : ! Intermediates
360 : integer :: i, j, counter
361 :
362 : ! Ensures that anything is less than initial i_min
363 : ! and anything is morre than initial i_max.
364 0 : do j=0,nz
365 0 : i_min(j) = nz+1
366 0 : i_max(j) = -1
367 : end do
368 :
369 : ! Find correct answers
370 0 : do counter=1,2*nz
371 0 : j = ranges(counter,1)
372 0 : i = ranges(counter,2)
373 0 : if (i < i_min(j)) i_min(j) = i
374 0 : if (i > i_max(j)) i_max(j) = i
375 : end do
376 :
377 0 : end subroutine find_i_ranges
378 :
379 : ! Computes the fraction of the mass which ends in cell j that passed through cell i
380 : ! for all (i,j) such that mesh_intersects(i,j) /= 0. These are the (i,j) for which
381 : ! the pass_fraction is not required to be 0 or 1, so we only need to precompute them.
382 : ! To do this we note that these (i,j) form a monotonic path from (0,0) to (nz,nz),
383 : ! and that within a row or column the path must be contiguous because cells represent
384 : ! contiguous mass intervals. This allows us to just sum along j in the direction the
385 : ! material flows to get the cumulative fraction of each cell which begins to one side
386 : ! of each face. Note that for cells where the flow converges (i.e. mass_flux(j) > 0,
387 : ! mass_flux(j+1) < 0) we perform a cumulative sum in both directions, meeting at
388 : ! cell j, at which point the two sums are added together.
389 0 : subroutine prepare_pass_fraction(nz, delta_m, dm, mesh_intersects, ranges, i_min, i_max, pf)
390 : ! Inputs
391 : integer :: nz
392 : real(qp) :: delta_m
393 : real(qp), dimension(:) :: dm, mesh_intersects
394 : integer, dimension(:,:) :: ranges
395 : integer, dimension(0:nz) :: i_min, i_max
396 : type(non_rect_array), dimension(0:nz) :: pf
397 :
398 : ! Intermediates
399 : integer :: i, j, l, counter
400 :
401 : ! Outputs
402 0 : call find_i_ranges(nz, ranges, i_min, i_max)
403 :
404 0 : do j=0,nz
405 0 : allocate(pf(j)%arr(i_min(j):i_max(j)))
406 0 : pf(j)%arr = 0
407 : end do
408 :
409 : ! Store mesh_intersects in pf
410 : counter = 1
411 0 : do j=0,nz
412 0 : do i=i_min(j),i_max(j)
413 0 : pf(j) % arr(i) = mesh_intersects(counter)
414 0 : counter = counter + 1
415 : end do
416 : end do
417 :
418 : ! Cumulatively sum mesh_intersects
419 : ! from i_min -> min(j, i_max) and from i_max -> max(j, i-min).
420 0 : do j=0,nz
421 0 : do i=i_min(j)+1,min(j,i_max(j))
422 0 : pf(j)%arr(i) = pf(j)%arr(i) + pf(j)%arr(i - 1)
423 : end do
424 0 : do i=i_max(j)-1,max(j,i_min(j)),-1
425 0 : pf(j)%arr(i) = pf(j)%arr(i) + pf(j)%arr(i + 1)
426 : end do
427 : end do
428 :
429 : ! Normalize by cell mass.
430 : ! Note that j == 0 only has a non-zero pass fraction when delta_m < 0,
431 : ! in which case it sums over l to -delta_m.
432 : ! So in all cases we may normalize the j == 0 case by -delta_m.
433 0 : do j=0,nz
434 0 : do l=i_min(j), i_max(j)
435 0 : if (j > 0) then
436 0 : pf(j)%arr(l) = pf(j)%arr(l) / dm(max(1,j)) ! bp: get rid of bogus compiler warning
437 : else
438 0 : pf(j)%arr(l) = -pf(j)%arr(l) / delta_m
439 : end if
440 : end do
441 : end do
442 :
443 0 : end subroutine prepare_pass_fraction
444 :
445 :
446 : ! Computes the fraction of mass in final cell j which was at some point during adjust_mass
447 : ! inside cell i.
448 0 : real(qp) function compute_pass_fraction(i, j, i_min, i_max, pf)
449 : ! Inputs
450 : integer :: i, j
451 : integer, dimension(0:) :: i_min, i_max
452 : type(non_rect_array), dimension(0:) :: pf
453 :
454 0 : if (i >= i_min(j) .and. i <= i_max(j)) then
455 0 : compute_pass_fraction = pf(j) % arr(i) ! In the pre-computed portion
456 0 : else if (i > max(j,i_max(j)) .or. i < min(j,i_min(j))) then
457 : compute_pass_fraction = 0 ! Outside of the flow
458 : else
459 0 : compute_pass_fraction = 1 ! In the flow but not pre-computed
460 : end if
461 :
462 0 : end function compute_pass_fraction
463 :
464 : ! Computes for each i the minimum and maximum j such that pass_fraction(i,j) > 0.
465 : ! The set of non-zero pass_fraction(i,:) is guaranteed to be contiguous so this
466 : ! is equivalent to enumerating all (i,j) with pass_fraction(i,j) > 0.
467 0 : subroutine find_j_ranges(nz, ranges, mass_flux, j_min, j_max)
468 : ! Inputs
469 : integer :: nz
470 : integer, dimension(:,:) :: ranges
471 : integer, dimension(:) :: j_min, j_max
472 : real(qp), dimension(:) :: mass_flux
473 :
474 : ! Intermediates
475 : integer :: i, j, counter
476 :
477 : ! Ensures that anything is less than initial i_min
478 : ! and anything is morre than initial i_max.
479 0 : do i=1,nz
480 0 : j_min(i) = nz+1
481 0 : j_max(i) = -1
482 : end do
483 :
484 : ! If the flow is downward then j_min = i (all lesser j have i > j).
485 : ! If the flow is downward then j_max is the greatest j such that (i,j) appears
486 : ! in ranges. All greater j have vanishing pass_fraction.
487 :
488 : ! If the flow is upward then j_max = i (all greater j have i < j).
489 : ! If the flow is upward then j_min is the least j such that (i,j) appears in ranges.
490 : ! All lesser j have vanishing pass_fraction.
491 :
492 : ! If material enters a cell from both faces then all material which touches
493 : ! the cell ends in the cell, so pass_fraction(i,:) is non-zero only for j == i.
494 : ! Hence j_min = j_max = i.
495 :
496 : ! If material leaves a cell from both faces then j_min is,
497 : ! the least j such that (i,j) appears in ranges, and j_max is
498 : ! the greatest j such that (i,j) appears in ranges.
499 :
500 : ! Note that for all of the above we may use inclusive inequalities
501 : ! because when one or more of the faces of a cell have zero mass flux
502 : ! these different cases agree.
503 :
504 : ! Set the easy ones
505 0 : do i=1,nz
506 0 : if (mass_flux(i) >= 0 .and. mass_flux(i+1) >= 0) then
507 0 : j_min(i) = i
508 : end if
509 :
510 0 : if (mass_flux(i) <= 0 .and. mass_flux(i+1) <= 0) then
511 0 : j_max(i) = i
512 : end if
513 :
514 0 : if (mass_flux(i) >= 0 .and. mass_flux(i+1) <= 0) then
515 0 : j_min(i) = i
516 0 : j_max(i) = i
517 : end if
518 :
519 : end do
520 :
521 : ! Now go through ranges
522 0 : do counter=1,2*nz
523 0 : j = ranges(counter,1)
524 0 : i = ranges(counter,2)
525 :
526 0 : if (i > 0) then
527 0 : if (mass_flux(i) >= 0 .and. mass_flux(i+1) >= 0) then
528 0 : if (j > j_max(i)) j_max(i) = j
529 : end if
530 :
531 0 : if (mass_flux(i) <= 0 .and. mass_flux(i+1) <= 0) then
532 0 : if (j < j_min(i)) j_min(i) = j
533 : end if
534 :
535 0 : if (mass_flux(i) <= 0 .and. mass_flux(i+1) >= 0) then
536 0 : if (j < j_min(i)) j_min(i) = j
537 0 : if (j > j_max(i)) j_max(i) = j
538 : end if
539 : end if
540 : end do
541 :
542 0 : end subroutine find_j_ranges
543 :
544 : ! Integrates conserved quantities given in vals_old from the mesh specified by dm_old onto the mesh given by dm_new.
545 : ! vals_new is set equal to \int vals_old dm / dm_new, where the integral runs over the mass range of the new cell.
546 : ! vals_outside is the value to use in the integral when the new cell runs outside the old mesh.
547 0 : subroutine integrate_conserved(vals_new, vals_old, vals_outside, dm_new, dm_old, nz, nvars)
548 :
549 : ! Inputs
550 : integer, intent(in) :: nz, nvars
551 : real(dp), intent(in), dimension(:,:) :: vals_old ! (cell index, variable index)
552 : real(dp), intent(in), dimension(:) :: vals_outside ! (variable index)
553 : real(qp), intent(in), dimension(:) :: dm_new, dm_old
554 :
555 : ! Intermediates
556 : integer :: i, j, k, l
557 0 : real(qp) :: mesh_intersects(2*nz)
558 0 : integer :: ranges(2*nz,2)
559 :
560 : ! Output
561 : real(dp), intent(out), dimension(:,:) :: vals_new
562 :
563 0 : ranges = 0
564 0 : mesh_intersects = 0d0
565 :
566 : ! Intersect meshes
567 0 : call make_compressed_intersect(dm_new, dm_old, nz, mesh_intersects, ranges)
568 :
569 0 : vals_new = 0d0
570 0 : do i=1,2*nz
571 0 : j = ranges(i,1)
572 0 : k = ranges(i,2)
573 :
574 0 : if (j == 0) cycle
575 :
576 0 : do l=1,nvars
577 0 : if (k > 0) then
578 0 : vals_new(j,l) = vals_new(j,l) + mesh_intersects(i) * vals_old(k,l)
579 : else
580 0 : vals_new(j,l) = vals_new(j,l) + mesh_intersects(i) * vals_outside(l)
581 : end if
582 : end do
583 :
584 : end do
585 :
586 0 : do j=1,nz
587 0 : do l=1,nvars
588 0 : vals_new(j,l) = vals_new(j,l) / dm_new(j)
589 : end do
590 : end do
591 :
592 0 : end subroutine integrate_conserved
593 :
594 0 : end module mass_utils
|
