PQM_functions.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Piecewise quartic reconstruction functions
6module pqm_functions
7
9
10implicit none ; private
11
13
14contains
15
16!> Reconstruction by quartic polynomials within each cell.
17!!
18!! It is assumed that the dimension of 'u' is equal to the number of cells
19!! defining 'grid' and 'ppoly'. No consistency check is performed.
20subroutine pqm_reconstruction( N, h, u, edge_values, edge_slopes, ppoly_coef, h_neglect, answer_date )
21 integer, intent(in) :: n !< Number of cells
22 real, dimension(:), intent(in) :: h !< cell widths (size N) [H]
23 real, dimension(:), intent(in) :: u !< cell averages (size N) [A]
24 real, dimension(:,:), intent(inout) :: edge_values !< Edge value of polynomial [A]
25 real, dimension(:,:), intent(inout) :: edge_slopes !< Edge slope of polynomial [A H-1]
26 real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial, mainly [A]
27 real, intent(in) :: h_neglect !< A negligibly small width for
28 !! the purpose of cell reconstructions [H]
29 integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use
30
31 ! Local variables
32 integer :: k ! loop index
33 real :: h_c ! cell width [H]
34 real :: u0_l, u0_r ! edge values (left and right) [A]
35 real :: u1_l, u1_r ! edge slopes (left and right) [A H-1]
36 real :: a, b, c, d, e ! quartic fit coefficients [A]
37
38 ! PQM limiter
39 call pqm_limiter( n, h, u, edge_values, edge_slopes, h_neglect, answer_date=answer_date )
40
41 ! Loop on cells to construct the cubic within each cell
42 do k = 1,n
43
44 u0_l = edge_values(k,1)
45 u0_r = edge_values(k,2)
46
47 u1_l = edge_slopes(k,1)
48 u1_r = edge_slopes(k,2)
49
50 h_c = h(k)
51
52 a = u0_l
53 b = h_c * u1_l
54 c = 30.0 * u(k) - 12.0*u0_r - 18.0*u0_l + 1.5*h_c*(u1_r - 3.0*u1_l)
55 d = -60.0 * u(k) + h_c *(6.0*u1_l - 4.0*u1_r) + 28.0*u0_r + 32.0*u0_l
56 e = 30.0 * u(k) + 2.5*h_c*(u1_r - u1_l) - 15.0*(u0_l + u0_r)
57
58 ! Store coefficients
59 ppoly_coef(k,1) = a
60 ppoly_coef(k,2) = b
61 ppoly_coef(k,3) = c
62 ppoly_coef(k,4) = d
63 ppoly_coef(k,5) = e
64
65 enddo ! end loop on cells
66
67end subroutine pqm_reconstruction
68
69!> Limit the piecewise quartic method reconstruction
70!!
71!! Standard PQM limiter (White & Adcroft, JCP 2008).
72!!
73!! It is assumed that the dimension of 'u' is equal to the number of cells
74!! defining 'grid' and 'ppoly'. No consistency check is performed.
75subroutine pqm_limiter( N, h, u, edge_values, edge_slopes, h_neglect, answer_date )
76 integer, intent(in) :: N !< Number of cells
77 real, dimension(:), intent(in) :: h !< cell widths (size N) [H]
78 real, dimension(:), intent(in) :: u !< cell average properties (size N) [A]
79 real, dimension(:,:), intent(inout) :: edge_values !< Potentially modified edge values [A]
80 real, dimension(:,:), intent(inout) :: edge_slopes !< Potentially modified edge slopes [A H-1]
81 real, intent(in) :: h_neglect !< A negligibly small width for
82 !! the purpose of cell reconstructions [H]
83 integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use
84
85 ! Local variables
86 integer :: k ! loop index
87 integer :: inflexion_l
88 integer :: inflexion_r
89 real :: u0_l, u0_r ! edge values [A]
90 real :: u1_l, u1_r ! edge slopes [A H-1]
91 real :: u_l, u_c, u_r ! left, center and right cell averages [A]
92 real :: h_l, h_c, h_r ! left, center and right cell widths [H]
93 real :: sigma_l, sigma_c, sigma_r ! left, center and right van Leer slopes [A H-1]
94 real :: slope ! retained PLM slope [A H-1]
95 real :: a, b, c, d, e ! quartic fit coefficients [A]
96 real :: alpha1, alpha2, alpha3 ! Normalized second derivative coefficients [A]
97 real :: rho ! A temporary expression [A2]
98 real :: sqrt_rho ! The square root of rho [A]
99 real :: gradient1, gradient2 ! Normalized gradients [A]
100 real :: x1, x2 ! Fractional inflection point positions in a cell [nondim]
101
102 ! Bound edge values
103 call bound_edge_values( n, h, u, edge_values, h_neglect, answer_date=answer_date )
104
105 ! Make discontinuous edge values monotonic (thru averaging)
106 call check_discontinuous_edge_values( n, u, edge_values )
107
108 ! Loop on interior cells to apply the PQM limiter
109 do k = 2,n-1
110
111 !if ( h(k) < 1.0 ) cycle
112
113 inflexion_l = 0
114 inflexion_r = 0
115
116 ! Get edge values, edge slopes and cell width
117 u0_l = edge_values(k,1)
118 u0_r = edge_values(k,2)
119 u1_l = edge_slopes(k,1)
120 u1_r = edge_slopes(k,2)
121
122 ! Get cell widths and cell averages (boundary cells are assumed to
123 ! be local extrema for the sake of slopes)
124 h_l = h(k-1)
125 h_c = h(k)
126 h_r = h(k+1)
127 u_l = u(k-1)
128 u_c = u(k)
129 u_r = u(k+1)
130
131 ! Compute limited slope
132 sigma_l = 2.0 * ( u_c - u_l ) / ( h_c + h_neglect )
133 sigma_c = 2.0 * ( u_r - u_l ) / ( h_l + 2.0*h_c + h_r + h_neglect )
134 sigma_r = 2.0 * ( u_r - u_c ) / ( h_c + h_neglect )
135
136 if ( (sigma_l * sigma_r) > 0.0 ) then
137 slope = sign( min(abs(sigma_l),abs(sigma_c),abs(sigma_r)), sigma_c )
138 else
139 slope = 0.0
140 endif
141
142 ! If one of the slopes has the wrong sign compared with the
143 ! limited PLM slope, it is set equal to the limited PLM slope
144 if ( u1_l*slope <= 0.0 ) u1_l = slope
145 if ( u1_r*slope <= 0.0 ) u1_r = slope
146
147 ! Local extremum --> flatten
148 if ( (u0_r - u_c) * (u_c - u0_l) <= 0.0) then
149 u0_l = u_c
150 u0_r = u_c
151 u1_l = 0.0
152 u1_r = 0.0
153 inflexion_l = -1
154 inflexion_r = -1
155 endif
156
157 ! Edge values are bounded and averaged when discontinuous and not
158 ! monotonic, edge slopes are consistent and the cell is not an extremum.
159 ! We now need to check and enforce the monotonicity of the quartic within
160 ! the cell
161 if ( (inflexion_l == 0) .AND. (inflexion_r == 0) ) then
162
163 a = u0_l
164 b = h_c * u1_l
165 c = 30.0 * u(k) - 12.0*u0_r - 18.0*u0_l + 1.5*h_c*(u1_r - 3.0*u1_l)
166 d = -60.0 * u(k) + h_c *(6.0*u1_l - 4.0*u1_r) + 28.0*u0_r + 32.0*u0_l
167 e = 30.0 * u(k) + 2.5*h_c*(u1_r - u1_l) - 15.0*(u0_l + u0_r)
168
169 ! Determine the coefficients of the second derivative
170 ! alpha1 xi^2 + alpha2 xi + alpha3
171 alpha1 = 6*e
172 alpha2 = 3*d
173 alpha3 = c
174
175 rho = alpha2 * alpha2 - 4.0 * alpha1 * alpha3
176
177 ! Check whether inflexion points exist
178 if (( alpha1 /= 0.0 ) .and. ( rho >= 0.0 )) then
179
180 sqrt_rho = sqrt( rho )
181
182 x1 = 0.5 * ( - alpha2 - sqrt_rho ) / alpha1
183 x2 = 0.5 * ( - alpha2 + sqrt_rho ) / alpha1
184
185 ! Check whether both inflexion points lie in [0,1]
186 if ( (x1 >= 0.0) .AND. (x1 <= 1.0) .AND. &
187 (x2 >= 0.0) .AND. (x2 <= 1.0) ) then
188
189 gradient1 = 4.0 * e * (x1**3) + 3.0 * d * (x1**2) + 2.0 * c * x1 + b
190 gradient2 = 4.0 * e * (x2**3) + 3.0 * d * (x2**2) + 2.0 * c * x2 + b
191
192 ! Check whether one of the gradients is inconsistent
193 if ( (gradient1 * slope < 0.0) .OR. &
194 (gradient2 * slope < 0.0) ) then
195 ! Decide where to collapse inflexion points
196 ! (depends on one-sided slopes)
197 if ( abs(sigma_l) < abs(sigma_r) ) then
198 inflexion_l = 1
199 else
200 inflexion_r = 1
201 endif
202 endif
203
204 ! If both x1 and x2 do not lie in [0,1], check whether
205 ! only x1 lies in [0,1]
206 elseif ( (x1 >= 0.0) .AND. (x1 <= 1.0) ) then
207
208 gradient1 = 4.0 * e * (x1**3) + 3.0 * d * (x1**2) + 2.0 * c * x1 + b
209
210 ! Check whether the gradient is inconsistent
211 if ( gradient1 * slope < 0.0 ) then
212 ! Decide where to collapse inflexion points
213 ! (depends on one-sided slopes)
214 if ( abs(sigma_l) < abs(sigma_r) ) then
215 inflexion_l = 1
216 else
217 inflexion_r = 1
218 endif
219 endif
220
221 ! If x1 does not lie in [0,1], check whether x2 lies in [0,1]
222 elseif ( (x2 >= 0.0) .AND. (x2 <= 1.0) ) then
223
224 gradient2 = 4.0 * e * (x2**3) + 3.0 * d * (x2**2) + 2.0 * c * x2 + b
225
226 ! Check whether the gradient is inconsistent
227 if ( gradient2 * slope < 0.0 ) then
228 ! Decide where to collapse inflexion points
229 ! (depends on one-sided slopes)
230 if ( abs(sigma_l) < abs(sigma_r) ) then
231 inflexion_l = 1
232 else
233 inflexion_r = 1
234 endif
235 endif
236
237 endif ! end checking where the inflexion points lie
238
239 endif ! end checking if alpha1 != 0 AND rho >= 0
240
241 ! If alpha1 is zero, the second derivative of the quartic reduces
242 ! to a straight line
243 if (( alpha1 == 0.0 ) .and. ( alpha2 /= 0.0 )) then
244
245 x1 = - alpha3 / alpha2
246 if ( (x1 >= 0.0) .AND. (x1 <= 1.0) ) then
247
248 gradient1 = 4.0 * e * (x1**3) + 3.0 * d * (x1**2) + 2.0 * c * x1 + b
249
250 ! Check whether the gradient is inconsistent
251 if ( gradient1 * slope < 0.0 ) then
252 ! Decide where to collapse inflexion points
253 ! (depends on one-sided slopes)
254 if ( abs(sigma_l) < abs(sigma_r) ) then
255 inflexion_l = 1
256 else
257 inflexion_r = 1
258 endif
259 endif ! check slope consistency
260
261 endif
262
263 endif ! end check whether we can find the root of the straight line
264
265 endif ! end checking whether to shift inflexion points
266
267 ! At this point, we know onto which edge to shift inflexion points
268 if ( inflexion_l == 1 ) then
269
270 ! We modify the edge slopes so that both inflexion points
271 ! collapse onto the left edge
272 u1_l = ( 10.0 * u_c - 2.0 * u0_r - 8.0 * u0_l ) / (3.0*h_c + h_neglect )
273 u1_r = ( -10.0 * u_c + 6.0 * u0_r + 4.0 * u0_l ) / ( h_c + h_neglect )
274
275 ! One of the modified slopes might be inconsistent. When that happens,
276 ! the inconsistent slope is set equal to zero and the opposite edge value
277 ! and edge slope are modified in compliance with the fact that both
278 ! inflexion points must still be located on the left edge
279 if ( u1_l * slope < 0.0 ) then
280
281 u1_l = 0.0
282 u0_r = 5.0 * u_c - 4.0 * u0_l
283 u1_r = 20.0 * (u_c - u0_l) / ( h_c + h_neglect )
284
285 elseif ( u1_r * slope < 0.0 ) then
286
287 u1_r = 0.0
288 u0_l = (5.0*u_c - 3.0*u0_r) / 2.0
289 u1_l = 10.0 * (-u_c + u0_r) / (3.0 * h_c + h_neglect)
290
291 endif
292
293 elseif ( inflexion_r == 1 ) then
294
295 ! We modify the edge slopes so that both inflexion points
296 ! collapse onto the right edge
297 u1_r = ( -10.0 * u_c + 8.0 * u0_r + 2.0 * u0_l ) / (3.0 * h_c + h_neglect)
298 u1_l = ( 10.0 * u_c - 4.0 * u0_r - 6.0 * u0_l ) / (h_c + h_neglect)
299
300 ! One of the modified slopes might be inconsistent. When that happens,
301 ! the inconsistent slope is set equal to zero and the opposite edge value
302 ! and edge slope are modified in compliance with the fact that both
303 ! inflexion points must still be located on the right edge
304 if ( u1_l * slope < 0.0 ) then
305
306 u1_l = 0.0
307 u0_r = ( 5.0 * u_c - 3.0 * u0_l ) / 2.0
308 u1_r = 10.0 * (u_c - u0_l) / (3.0 * h_c + h_neglect)
309
310 elseif ( u1_r * slope < 0.0 ) then
311
312 u1_r = 0.0
313 u0_l = 5.0 * u_c - 4.0 * u0_r
314 u1_l = 20.0 * ( -u_c + u0_r ) / (h_c + h_neglect)
315
316 endif
317
318 endif ! clause to check where to collapse inflexion points
319
320 ! Save edge values and edge slopes for reconstruction
321 edge_values(k,1) = u0_l
322 edge_values(k,2) = u0_r
323 edge_slopes(k,1) = u1_l
324 edge_slopes(k,2) = u1_r
325
326 enddo ! end loop on interior cells
327
328 ! Constant reconstruction within boundary cells
329 edge_values(1,:) = u(1)
330 edge_slopes(1,:) = 0.0
331
332 edge_values(n,:) = u(n)
333 edge_slopes(n,:) = 0.0
334
335end subroutine pqm_limiter
336
337!> Reconstruction by parabolas within boundary cells.
338!!
339!! The following explanations apply to the left boundary cell. The same
340!! reasoning holds for the right boundary cell.
341!!
342!! A parabola needs to be built in the cell and requires three degrees of
343!! freedom, which are the right edge value and slope and the cell average.
344!! The right edge values and slopes are taken to be that of the neighboring
345!! cell (i.e., the left edge value and slope of the neighboring cell).
346!! The resulting parabola is not necessarily monotonic and the traditional
347!! PPM limiter is used to modify one of the edge values in order to yield
348!! a monotonic parabola.
349!!
350!! It is assumed that the size of the array 'u' is equal to the number of cells
351!! defining 'grid' and 'ppoly'. No consistency check is performed here.
352subroutine pqm_boundary_extrapolation( N, h, u, edge_values, ppoly_coef )
353 integer, intent(in) :: n !< Number of cells
354 real, dimension(:), intent(in) :: h !< cell widths (size N) [H]
355 real, dimension(:), intent(in) :: u !< cell averages (size N) [A]
356 real, dimension(:,:), intent(inout) :: edge_values !< Edge value of polynomial [A]
357 real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial, mainly [A]
358 ! Local variables
359 integer :: i0, i1
360 real :: u0, u1 ! Successive cell averages [A]
361 real :: h0, h1 ! Successive cell thicknesses [H]
362 real :: a, b, c, d, e ! quartic fit coefficients [A]
363 real :: u0_l, u0_r ! Edge values [A]
364 real :: u1_l, u1_r ! Edge slopes [A H-1]
365 real :: slope ! The integrated slope across the cell [A]
366 real :: exp1, exp2 ! Two temporary expressions [A2]
367
368 ! ----- Left boundary -----
369 i0 = 1
370 i1 = 2
371 h0 = h(i0)
372 h1 = h(i1)
373 u0 = u(i0)
374 u1 = u(i1)
375
376 ! Compute the left edge slope in neighboring cell and express it in
377 ! the global coordinate system
378 b = ppoly_coef(i1,2)
379 u1_r = b *(h0/h1) ! derivative evaluated at xi = 0.0,
380 ! expressed w.r.t. xi (local coord. system)
381
382 ! Limit the right slope by the PLM limited slope
383 slope = 2.0 * ( u1 - u0 )
384 if ( abs(u1_r) > abs(slope) ) then
385 u1_r = slope
386 endif
387
388 ! The right edge value in the boundary cell is taken to be the left
389 ! edge value in the neighboring cell
390 u0_r = edge_values(i1,1)
391
392 ! Given the right edge value and slope, we determine the left
393 ! edge value and slope by computing the parabola as determined by
394 ! the right edge value and slope and the boundary cell average
395 u0_l = 3.0 * u0 + 0.5 * u1_r - 2.0 * u0_r
396
397 ! Apply the traditional PPM limiter
398 exp1 = (u0_r - u0_l) * (u0 - 0.5*(u0_l+u0_r))
399 exp2 = (u0_r - u0_l) * (u0_r - u0_l) / 6.0
400
401 if ( exp1 > exp2 ) then
402 u0_l = 3.0 * u0 - 2.0 * u0_r
403 endif
404
405 if ( exp1 < -exp2 ) then
406 u0_r = 3.0 * u0 - 2.0 * u0_l
407 endif
408
409 edge_values(i0,1) = u0_l
410 edge_values(i0,2) = u0_r
411
412 a = u0_l
413 b = 6.0 * u0 - 4.0 * u0_l - 2.0 * u0_r
414 c = 3.0 * ( u0_r + u0_l - 2.0 * u0 )
415
416 ! The quartic is reduced to a parabola in the boundary cell
417 ppoly_coef(i0,1) = a
418 ppoly_coef(i0,2) = b
419 ppoly_coef(i0,3) = c
420 ppoly_coef(i0,4) = 0.0
421 ppoly_coef(i0,5) = 0.0
422
423 ! ----- Right boundary -----
424 i0 = n-1
425 i1 = n
426 h0 = h(i0)
427 h1 = h(i1)
428 u0 = u(i0)
429 u1 = u(i1)
430
431 ! Compute the right edge slope in neighboring cell and express it in
432 ! the global coordinate system
433 b = ppoly_coef(i0,2)
434 c = ppoly_coef(i0,3)
435 d = ppoly_coef(i0,4)
436 e = ppoly_coef(i0,5)
437 u1_l = (b + 2*c + 3*d + 4*e) ! derivative evaluated at xi = 1.0
438 u1_l = u1_l * (h1/h0)
439
440 ! Limit the left slope by the PLM limited slope
441 slope = 2.0 * ( u1 - u0 )
442 if ( abs(u1_l) > abs(slope) ) then
443 u1_l = slope
444 endif
445
446 ! The left edge value in the boundary cell is taken to be the right
447 ! edge value in the neighboring cell
448 u0_l = edge_values(i0,2)
449
450 ! Given the left edge value and slope, we determine the right
451 ! edge value and slope by computing the parabola as determined by
452 ! the left edge value and slope and the boundary cell average
453 u0_r = 3.0 * u1 - 0.5 * u1_l - 2.0 * u0_l
454
455 ! Apply the traditional PPM limiter
456 exp1 = (u0_r - u0_l) * (u1 - 0.5*(u0_l+u0_r))
457 exp2 = (u0_r - u0_l) * (u0_r - u0_l) / 6.0
458
459 if ( exp1 > exp2 ) then
460 u0_l = 3.0 * u1 - 2.0 * u0_r
461 endif
462
463 if ( exp1 < -exp2 ) then
464 u0_r = 3.0 * u1 - 2.0 * u0_l
465 endif
466
467 edge_values(i1,1) = u0_l
468 edge_values(i1,2) = u0_r
469
470 a = u0_l
471 b = 6.0 * u1 - 4.0 * u0_l - 2.0 * u0_r
472 c = 3.0 * ( u0_r + u0_l - 2.0 * u1 )
473
474 ! The quartic is reduced to a parabola in the boundary cell
475 ppoly_coef(i1,1) = a
476 ppoly_coef(i1,2) = b
477 ppoly_coef(i1,3) = c
478 ppoly_coef(i1,4) = 0.0
479 ppoly_coef(i1,5) = 0.0
480
481end subroutine pqm_boundary_extrapolation
482
483
484!> Reconstruction by parabolas within boundary cells.
485!!
486!! The following explanations apply to the left boundary cell. The same
487!! reasoning holds for the right boundary cell.
488!!
489!! A parabola needs to be built in the cell and requires three degrees of
490!! freedom, which are the right edge value and slope and the cell average.
491!! The right edge values and slopes are taken to be that of the neighboring
492!! cell (i.e., the left edge value and slope of the neighboring cell).
493!! The resulting parabola is not necessarily monotonic and the traditional
494!! PPM limiter is used to modify one of the edge values in order to yield
495!! a monotonic parabola.
496!!
497!! It is assumed that the size of the array 'u' is equal to the number of cells
498!! defining 'grid' and 'ppoly'. No consistency check is performed here.
499subroutine pqm_boundary_extrapolation_v1( N, h, u, edge_values, edge_slopes, ppoly_coef, h_neglect )
500 integer, intent(in) :: n !< Number of cells
501 real, dimension(:), intent(in) :: h !< cell widths (size N) [H]
502 real, dimension(:), intent(in) :: u !< cell averages (size N) [A]
503 real, dimension(:,:), intent(inout) :: edge_values !< Edge value of polynomial [A]
504 real, dimension(:,:), intent(inout) :: edge_slopes !< Edge slope of polynomial [A H-1]
505 real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial, mainly [A]
506 real, intent(in) :: h_neglect !< A negligibly small width for
507 !! the purpose of cell reconstructions [H]
508 ! Local variables
509 integer :: i0, i1
510 integer :: inflexion_l
511 integer :: inflexion_r
512 real :: u0, u1, um ! Successive cell averages [A]
513 real :: h0, h1 ! Successive cell thicknesses [H]
514 real :: a, b, c, d, e ! quartic fit coefficients [A]
515 real :: ar, br ! Temporary variables in [A]
516 real :: beta ! A rational function coefficient [nondim]
517 real :: u0_l, u0_r ! Edge values [A]
518 real :: u1_l, u1_r ! Edge slopes [A H-1]
519 real :: u_plm ! The integrated piecewise linear method slope [A]
520 real :: slope ! The integrated slope across the cell [A]
521 real :: alpha1, alpha2, alpha3 ! Normalized second derivative coefficients [A]
522 real :: rho ! A temporary expression [A2]
523 real :: sqrt_rho ! The square root of rho [A]
524 real :: gradient1, gradient2 ! Normalized gradients [A]
525 real :: x1, x2 ! Fractional inflection point positions in a cell [nondim]
526
527 ! ----- Left boundary (TOP) -----
528 i0 = 1
529 i1 = 2
530 h0 = h(i0)
531 h1 = h(i1)
532 u0 = u(i0)
533 u1 = u(i1)
534 um = u0
535
536 ! Compute real slope and express it w.r.t. local coordinate system
537 ! within boundary cell
538 slope = 2.0 * ( u1 - u0 ) / ( ( h0 + h1 ) + h_neglect )
539 slope = slope * h0
540
541 ! The right edge value and slope of the boundary cell are taken to be the
542 ! left edge value and slope of the adjacent cell
543 a = ppoly_coef(i1,1)
544 b = ppoly_coef(i1,2)
545
546 u0_r = a ! edge value
547 u1_r = b / (h1 + h_neglect) ! edge slope (w.r.t. global coord.)
548
549 ! Compute coefficient for rational function based on mean and right
550 ! edge value and slope
551 if (u1_r /= 0.) then ! HACK by AJA
552 beta = 2.0 * ( u0_r - um ) / ( (h0 + h_neglect)*u1_r) - 1.0
553 else
554 beta = 0.
555 endif ! HACK by AJA
556 br = u0_r + beta*u0_r - um
557 ar = um + beta*um - br
558
559 ! Left edge value estimate based on rational function
560 u0_l = ar
561
562 ! Edge value estimate based on PLM
563 u_plm = um - 0.5 * slope
564
565 ! Check whether the left edge value is bounded by the mean and
566 ! the PLM edge value. If so, keep it and compute left edge slope
567 ! based on the rational function. If not, keep the PLM edge value and
568 ! compute corresponding slope.
569 if ( abs(um-u0_l) < abs(um-u_plm) ) then
570 u1_l = 2.0 * ( br - ar*beta)
571 u1_l = u1_l / (h0 + h_neglect)
572 else
573 u0_l = u_plm
574 u1_l = slope / (h0 + h_neglect)
575 endif
576
577 ! Monotonize quartic
578 inflexion_l = 0
579
580 a = u0_l
581 b = h0 * u1_l
582 c = 30.0 * um - 12.0*u0_r - 18.0*u0_l + 1.5*h0*(u1_r - 3.0*u1_l)
583 d = -60.0 * um + h0 *(6.0*u1_l - 4.0*u1_r) + 28.0*u0_r + 32.0*u0_l
584 e = 30.0 * um + 2.5*h0*(u1_r - u1_l) - 15.0*(u0_l + u0_r)
585
586 alpha1 = 6*e
587 alpha2 = 3*d
588 alpha3 = c
589
590 rho = alpha2 * alpha2 - 4.0 * alpha1 * alpha3
591
592 ! Check whether inflexion points exist. If so, transform the quartic
593 ! so that both inflexion points coalesce on the left edge.
594 if (( alpha1 /= 0.0 ) .and. ( rho >= 0.0 )) then
595
596 sqrt_rho = sqrt( rho )
597
598 x1 = 0.5 * ( - alpha2 - sqrt_rho ) / alpha1
599 if ( (x1 > 0.0) .and. (x1 < 1.0) ) then
600 gradient1 = 4.0 * e * (x1**3) + 3.0 * d * (x1**2) + 2.0 * c * x1 + b
601 if ( gradient1 * slope < 0.0 ) then
602 inflexion_l = 1
603 endif
604 endif
605
606 x2 = 0.5 * ( - alpha2 + sqrt_rho ) / alpha1
607 if ( (x2 > 0.0) .and. (x2 < 1.0) ) then
608 gradient2 = 4.0 * e * (x2**3) + 3.0 * d * (x2**2) + 2.0 * c * x2 + b
609 if ( gradient2 * slope < 0.0 ) then
610 inflexion_l = 1
611 endif
612 endif
613
614 endif
615
616 if (( alpha1 == 0.0 ) .and. ( alpha2 /= 0.0 )) then
617
618 x1 = - alpha3 / alpha2
619 if ( (x1 >= 0.0) .and. (x1 <= 1.0) ) then
620 gradient1 = 3.0 * d * (x1**2) + 2.0 * c * x1 + b
621 if ( gradient1 * slope < 0.0 ) then
622 inflexion_l = 1
623 endif
624 endif
625
626 endif
627
628 if ( inflexion_l == 1 ) then
629
630 ! We modify the edge slopes so that both inflexion points
631 ! collapse onto the left edge
632 u1_l = ( 10.0 * um - 2.0 * u0_r - 8.0 * u0_l ) / (3.0*h0 + h_neglect)
633 u1_r = ( -10.0 * um + 6.0 * u0_r + 4.0 * u0_l ) / (h0 + h_neglect)
634
635 ! One of the modified slopes might be inconsistent. When that happens,
636 ! the inconsistent slope is set equal to zero and the opposite edge value
637 ! and edge slope are modified in compliance with the fact that both
638 ! inflexion points must still be located on the left edge
639 if ( u1_l * slope < 0.0 ) then
640
641 u1_l = 0.0
642 u0_r = 5.0 * um - 4.0 * u0_l
643 u1_r = 20.0 * (um - u0_l) / ( h0 + h_neglect )
644
645 elseif ( u1_r * slope < 0.0 ) then
646
647 u1_r = 0.0
648 u0_l = (5.0*um - 3.0*u0_r) / 2.0
649 u1_l = 10.0 * (-um + u0_r) / (3.0 * h0 + h_neglect )
650
651 endif
652
653 endif
654
655 ! Store edge values, edge slopes and coefficients
656 edge_values(i0,1) = u0_l
657 edge_values(i0,2) = u0_r
658 edge_slopes(i0,1) = u1_l
659 edge_slopes(i0,2) = u1_r
660
661 a = u0_l
662 b = h0 * u1_l
663 c = 30.0 * um - 12.0*u0_r - 18.0*u0_l + 1.5*h0*(u1_r - 3.0*u1_l)
664 d = -60.0 * um + h0 *(6.0*u1_l - 4.0*u1_r) + 28.0*u0_r + 32.0*u0_l
665 e = 30.0 * um + 2.5*h0*(u1_r - u1_l) - 15.0*(u0_l + u0_r)
666
667 ! Store coefficients
668 ppoly_coef(i0,1) = a
669 ppoly_coef(i0,2) = b
670 ppoly_coef(i0,3) = c
671 ppoly_coef(i0,4) = d
672 ppoly_coef(i0,5) = e
673
674 ! ----- Right boundary (BOTTOM) -----
675 i0 = n-1
676 i1 = n
677 h0 = h(i0)
678 h1 = h(i1)
679 u0 = u(i0)
680 u1 = u(i1)
681 um = u1
682
683 ! Compute real slope and express it w.r.t. local coordinate system
684 ! within boundary cell
685 slope = 2.0 * ( u1 - u0 ) / ( h0 + h1 )
686 slope = slope * h1
687
688 ! The left edge value and slope of the boundary cell are taken to be the
689 ! right edge value and slope of the adjacent cell
690 a = ppoly_coef(i0,1)
691 b = ppoly_coef(i0,2)
692 c = ppoly_coef(i0,3)
693 d = ppoly_coef(i0,4)
694 e = ppoly_coef(i0,5)
695 u0_l = a + b + c + d + e ! edge value
696 u1_l = (b + 2*c + 3*d + 4*e) / h0 ! edge slope (w.r.t. global coord.)
697
698 ! Compute coefficient for rational function based on mean and left
699 ! edge value and slope
700 if (um-u0_l /= 0.) then ! HACK by AJA
701 beta = 0.5*h1*u1_l / (um-u0_l) - 1.0
702 else
703 beta = 0.
704 endif ! HACK by AJA
705 br = beta*um + um - u0_l
706 ar = u0_l
707
708 ! Right edge value estimate based on rational function
709 if (1+beta /= 0.) then ! HACK by AJA
710 u0_r = (ar + 2*br + beta*br ) / ((1+beta)*(1+beta))
711 else
712 u0_r = um + 0.5 * slope ! PLM
713 endif ! HACK by AJA
714
715 ! Right edge value estimate based on PLM
716 u_plm = um + 0.5 * slope
717
718 ! Check whether the right edge value is bounded by the mean and
719 ! the PLM edge value. If so, keep it and compute right edge slope
720 ! based on the rational function. If not, keep the PLM edge value and
721 ! compute corresponding slope.
722 if ( abs(um-u0_r) < abs(um-u_plm) ) then
723 u1_r = 2.0 * ( br - ar*beta ) / ( (1+beta)*(1+beta)*(1+beta) )
724 u1_r = u1_r / h1
725 else
726 u0_r = u_plm
727 u1_r = slope / h1
728 endif
729
730 ! Monotonize quartic
731 inflexion_r = 0
732
733 a = u0_l
734 b = h1 * u1_l
735 c = 30.0 * um - 12.0*u0_r - 18.0*u0_l + 1.5*h1*(u1_r - 3.0*u1_l)
736 d = -60.0 * um + h1*(6.0*u1_l - 4.0*u1_r) + 28.0*u0_r + 32.0*u0_l
737 e = 30.0 * um + 2.5*h1*(u1_r - u1_l) - 15.0*(u0_l + u0_r)
738
739 alpha1 = 6*e
740 alpha2 = 3*d
741 alpha3 = c
742
743 rho = alpha2 * alpha2 - 4.0 * alpha1 * alpha3
744
745 ! Check whether inflexion points exist. If so, transform the quartic
746 ! so that both inflexion points coalesce on the right edge.
747 if (( alpha1 /= 0.0 ) .and. ( rho >= 0.0 )) then
748
749 sqrt_rho = sqrt( rho )
750
751 x1 = 0.5 * ( - alpha2 - sqrt_rho ) / alpha1
752 if ( (x1 > 0.0) .and. (x1 < 1.0) ) then
753 gradient1 = 4.0 * e * (x1**3) + 3.0 * d * (x1**2) + 2.0 * c * x1 + b
754 if ( gradient1 * slope < 0.0 ) then
755 inflexion_r = 1
756 endif
757 endif
758
759 x2 = 0.5 * ( - alpha2 + sqrt_rho ) / alpha1
760 if ( (x2 > 0.0) .and. (x2 < 1.0) ) then
761 gradient2 = 4.0 * e * (x2**3) + 3.0 * d * (x2**2) + 2.0 * c * x2 + b
762 if ( gradient2 * slope < 0.0 ) then
763 inflexion_r = 1
764 endif
765 endif
766
767 endif
768
769 if (( alpha1 == 0.0 ) .and. ( alpha2 /= 0.0 )) then
770
771 x1 = - alpha3 / alpha2
772 if ( (x1 >= 0.0) .and. (x1 <= 1.0) ) then
773 gradient1 = 3.0 * d * (x1**2) + 2.0 * c * x1 + b
774 if ( gradient1 * slope < 0.0 ) then
775 inflexion_r = 1
776 endif
777 endif
778
779 endif
780
781 if ( inflexion_r == 1 ) then
782
783 ! We modify the edge slopes so that both inflexion points
784 ! collapse onto the right edge
785 u1_r = ( -10.0 * um + 8.0 * u0_r + 2.0 * u0_l ) / (3.0 * h1)
786 u1_l = ( 10.0 * um - 4.0 * u0_r - 6.0 * u0_l ) / h1
787
788 ! One of the modified slopes might be inconsistent. When that happens,
789 ! the inconsistent slope is set equal to zero and the opposite edge value
790 ! and edge slope are modified in compliance with the fact that both
791 ! inflexion points must still be located on the right edge
792 if ( u1_l * slope < 0.0 ) then
793
794 u1_l = 0.0
795 u0_r = ( 5.0 * um - 3.0 * u0_l ) / 2.0
796 u1_r = 10.0 * (um - u0_l) / (3.0 * h1)
797
798 elseif ( u1_r * slope < 0.0 ) then
799
800 u1_r = 0.0
801 u0_l = 5.0 * um - 4.0 * u0_r
802 u1_l = 20.0 * ( -um + u0_r ) / h1
803
804 endif
805
806 endif
807
808 ! Store edge values, edge slopes and coefficients
809 edge_values(i1,1) = u0_l
810 edge_values(i1,2) = u0_r
811 edge_slopes(i1,1) = u1_l
812 edge_slopes(i1,2) = u1_r
813
814 a = u0_l
815 b = h1 * u1_l
816 c = 30.0 * um - 12.0*u0_r - 18.0*u0_l + 1.5*h1*(u1_r - 3.0*u1_l)
817 d = -60.0 * um + h1 *(6.0*u1_l - 4.0*u1_r) + 28.0*u0_r + 32.0*u0_l
818 e = 30.0 * um + 2.5*h1*(u1_r - u1_l) - 15.0*(u0_l + u0_r)
819
820 ppoly_coef(i1,1) = a
821 ppoly_coef(i1,2) = b
822 ppoly_coef(i1,3) = c
823 ppoly_coef(i1,4) = d
824 ppoly_coef(i1,5) = e
825
827
828!> \namespace pqm_functions
829!!
830!! Date of creation: 2008.06.06
831!! L. White
832!!
833!! This module contains routines that handle one-dimensional finite volume
834!! reconstruction using the piecewise quartic method (PQM).
835
836end module pqm_functions