MOM6
P3M_functions.F90
Go to the documentation of this file.
1 !> Cubic interpolation functions
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
7 
8 implicit none ; private
9 
10 public p3m_interpolation
12 
13 real, parameter :: hneglect_dflt = 1.e-30 !< Default value of a negligible cell thickness
14 real, parameter :: hneglect_edge_dflt = 1.e-10 !< Default value of a negligible edge thickness
15 
16 contains
17 
18 !> Set up a piecewise cubic interpolation from cell averages and estimated
19 !! edge slopes and values
20 !!
21 !! Cubic interpolation between edges.
22 !!
23 !! The edge values and slopes MUST have been estimated prior to calling
24 !! this routine.
25 !!
26 !! It is assumed that the size of the array 'u' is equal to the number of cells
27 !! defining 'grid' and 'ppoly'. No consistency check is performed here.
28 subroutine p3m_interpolation( N, h, u, ppoly_E, ppoly_S, ppoly_coef, h_neglect )
29  integer, intent(in) :: n !< Number of cells
30  real, dimension(:), intent(in) :: h !< cell widths (size N) [H]
31  real, dimension(:), intent(in) :: u !< cell averages (size N) in arbitrary units [A]
32  real, dimension(:,:), intent(inout) :: ppoly_e !< Edge value of polynomial [A]
33  real, dimension(:,:), intent(inout) :: ppoly_s !< Edge slope of polynomial [A H-1].
34  real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial [A]
35  real, optional, intent(in) :: h_neglect !< A negligibly small width for the
36  !! purpose of cell reconstructions [H]
37 
38  ! Call the limiter for p3m, which takes care of everything from
39  ! computing the coefficients of the cubic to monotonizing it.
40  ! This routine could be called directly instead of having to call
41  ! 'P3M_interpolation' first but we do that to provide an homogeneous
42  ! interface.
43  call p3m_limiter( n, h, u, ppoly_e, ppoly_s, ppoly_coef, h_neglect )
44 
45 end subroutine p3m_interpolation
46 
47 !> Adust a piecewise cubic reconstruction with a limiter that adjusts the edge
48 !! values and slopes
49 !!
50 !! The p3m limiter operates as follows:
51 !!
52 !! 1. Edge values are bounded
53 !! 2. Discontinuous edge values are systematically averaged
54 !! 3. Loop on cells and do the following
55 !! a. Build cubic curve
56 !! b. Check if cubic curve is monotonic
57 !! c. If not, monotonize cubic curve and rebuild it
58 !!
59 !! Step 3 of the monotonization process leaves all edge values unchanged.
60 subroutine p3m_limiter( N, h, u, ppoly_E, ppoly_S, ppoly_coef, h_neglect )
61  integer, intent(in) :: N !< Number of cells
62  real, dimension(:), intent(in) :: h !< cell widths (size N) [H]
63  real, dimension(:), intent(in) :: u !< cell averages (size N) in arbitrary units [A]
64  real, dimension(:,:), intent(inout) :: ppoly_E !< Edge value of polynomial [A]
65  real, dimension(:,:), intent(inout) :: ppoly_S !< Edge slope of polynomial [A H-1]
66  real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial [A]
67  real, optional, intent(in) :: h_neglect !< A negligibly small width for
68  !! the purpose of cell reconstructions [H]
69  ! Local variables
70  integer :: k ! loop index
71  logical :: monotonic ! boolean indicating whether the cubic is monotonic
72  real :: u0_l, u0_r ! edge values [A]
73  real :: u1_l, u1_r ! edge slopes [A H-1]
74  real :: u_l, u_c, u_r ! left, center and right cell averages [A]
75  real :: h_l, h_c, h_r ! left, center and right cell widths [H]
76  real :: sigma_l, sigma_c, sigma_r ! left, center and right van Leer slopes [A H-1]
77  real :: slope ! retained PLM slope [A H-1]
78  real :: eps
79  real :: hNeglect ! A negligibly small thickness [H]
80 
81  hneglect = hneglect_dflt ; if (present(h_neglect)) hneglect = h_neglect
82 
83  eps = 1e-10
84 
85  ! 1. Bound edge values (boundary cells are assumed to be local extrema)
86  call bound_edge_values( n, h, u, ppoly_e, hneglect )
87 
88  ! 2. Systematically average discontinuous edge values
89  call average_discontinuous_edge_values( n, ppoly_e )
90 
91 
92  ! 3. Loop on cells and do the following
93  ! (a) Build cubic curve
94  ! (b) Check if cubic curve is monotonic
95  ! (c) If not, monotonize cubic curve and rebuild it
96  do k = 1,n
97 
98  ! Get edge values, edge slopes and cell width
99  u0_l = ppoly_e(k,1)
100  u0_r = ppoly_e(k,2)
101  u1_l = ppoly_s(k,1)
102  u1_r = ppoly_s(k,2)
103 
104  ! Get cell widths and cell averages (boundary cells are assumed to
105  ! be local extrema for the sake of slopes)
106  u_c = u(k)
107  h_c = h(k)
108 
109  if ( k == 1 ) then
110  h_l = h(k)
111  u_l = u(k)
112  else
113  h_l = h(k-1)
114  u_l = u(k-1)
115  endif
116 
117  if ( k == n ) then
118  h_r = h(k)
119  u_r = u(k)
120  else
121  h_r = h(k+1)
122  u_r = u(k+1)
123  endif
124 
125  ! Compute limited slope
126  sigma_l = 2.0 * ( u_c - u_l ) / ( h_c + hneglect )
127  sigma_c = 2.0 * ( u_r - u_l ) / ( h_l + 2.0*h_c + h_r + hneglect )
128  sigma_r = 2.0 * ( u_r - u_c ) / ( h_c + hneglect )
129 
130  if ( (sigma_l * sigma_r) > 0.0 ) then
131  slope = sign( min(abs(sigma_l),abs(sigma_c),abs(sigma_r)), sigma_c )
132  else
133  slope = 0.0
134  endif
135 
136  ! If the slopes are small, set them to zero to prevent asymmetric representation near extrema.
137  if ( abs(u1_l*h_c) < epsilon(u_c)*abs(u_c) ) u1_l = 0.0
138  if ( abs(u1_r*h_c) < epsilon(u_c)*abs(u_c) ) u1_r = 0.0
139 
140  ! The edge slopes are limited from above by the respective
141  ! one-sided slopes
142  if ( abs(u1_l) > abs(sigma_l) ) then
143  u1_l = sigma_l
144  endif
145 
146  if ( abs(u1_r) > abs(sigma_r) ) then
147  u1_r = sigma_r
148  endif
149 
150  ! Build cubic interpolant (compute the coefficients)
151  call build_cubic_interpolant( h, k, ppoly_e, ppoly_s, ppoly_coef )
152 
153  ! Check whether cubic is monotonic
154  monotonic = is_cubic_monotonic( ppoly_coef, k )
155 
156  ! If cubic is not monotonic, monotonize it by modifiying the
157  ! edge slopes, store the new edge slopes and recompute the
158  ! cubic coefficients
159  if ( .not.monotonic ) then
160  call monotonize_cubic( h_c, u0_l, u0_r, sigma_l, sigma_r, slope, u1_l, u1_r )
161  endif
162 
163  ! Store edge slopes
164  ppoly_s(k,1) = u1_l
165  ppoly_s(k,2) = u1_r
166 
167  ! Recompute coefficients of cubic
168  call build_cubic_interpolant( h, k, ppoly_e, ppoly_s, ppoly_coef )
169 
170  enddo ! loop on cells
171 
172 end subroutine p3m_limiter
173 
174 
175 !> Calculate the edge values and slopes at boundary cells as part of building a
176 !! piecewise cubic sub-grid scale profiles
177 !!
178 !! The following explanations apply to the left boundary cell. The same
179 !! reasoning holds for the right boundary cell.
180 !!
181 !! A cubic needs to be built in the cell and requires four degrees of freedom,
182 !! which are the edge values and slopes. The right edge values and slopes are
183 !! taken to be that of the neighboring cell (i.e., the left edge value and slope
184 !! of the neighboring cell). The left edge value and slope are determined by
185 !! computing the parabola based on the cell average and the right edge value
186 !! and slope. The resulting cubic is not necessarily monotonic and the slopes
187 !! are subsequently modified to yield a monotonic cubic.
188 subroutine p3m_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coef, &
189  h_neglect, h_neglect_edge )
190  integer, intent(in) :: n !< Number of cells
191  real, dimension(:), intent(in) :: h !< cell widths (size N) [H]
192  real, dimension(:), intent(in) :: u !< cell averages (size N) in arbitrary units [A]
193  real, dimension(:,:), intent(inout) :: ppoly_e !< Edge value of polynomial [A]
194  real, dimension(:,:), intent(inout) :: ppoly_s !< Edge slope of polynomial [A H-1]
195  real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial [A]
196  real, optional, intent(in) :: h_neglect !< A negligibly small width for the
197  !! purpose of cell reconstructions [H]
198  real, optional, intent(in) :: h_neglect_edge !< A negligibly small width
199  !! for the purpose of finding edge values [H]
200  ! Local variables
201  integer :: i0, i1
202  logical :: monotonic ! boolean indicating whether the cubic is monotonic
203  real :: u0, u1 ! Values of u in two adjacent cells [A]
204  real :: h0, h1 ! Values of h in two adjacent cells, plus a smal increment [H]
205  real :: b, c, d ! Temporary variables [A]
206  real :: u0_l, u0_r ! Left and right edge values [A]
207  real :: u1_l, u1_r ! Left and right edge slopes [A H-1]
208  real :: slope ! The cell center slope [A H-1]
209  real :: hneglect, hneglect_edge ! Negligibly small thickness [H]
210 
211  hneglect = hneglect_dflt ; if (present(h_neglect)) hneglect = h_neglect
212  hneglect_edge = hneglect_edge_dflt
213  if (present(h_neglect_edge)) hneglect_edge = h_neglect_edge
214 
215  ! ----- Left boundary -----
216  i0 = 1
217  i1 = 2
218  h0 = h(i0) + hneglect_edge
219  h1 = h(i1) + hneglect_edge
220  u0 = u(i0)
221  u1 = u(i1)
222 
223  ! Compute the left edge slope in neighboring cell and express it in
224  ! the global coordinate system
225  b = ppoly_coef(i1,2)
226  u1_r = b / h1 ! derivative evaluated at xi = 0.0, expressed w.r.t. x
227 
228  ! Limit the right slope by the PLM limited slope
229  slope = 2.0 * ( u1 - u0 ) / ( h0 + hneglect )
230  if ( abs(u1_r) > abs(slope) ) then
231  u1_r = slope
232  endif
233 
234  ! The right edge value in the boundary cell is taken to be the left
235  ! edge value in the neighboring cell
236  u0_r = ppoly_e(i1,1)
237 
238  ! Given the right edge value and slope, we determine the left
239  ! edge value and slope by computing the parabola as determined by
240  ! the right edge value and slope and the boundary cell average
241  u0_l = 3.0 * u0 + 0.5 * h0*u1_r - 2.0 * u0_r
242  u1_l = ( - 6.0 * u0 - 2.0 * h0*u1_r + 6.0 * u0_r) / ( h0 + hneglect )
243 
244  ! Check whether the edge values are monotonic. For example, if the left edge
245  ! value is larger than the right edge value while the slope is positive, the
246  ! edge values are inconsistent and we need to modify the left edge value
247  if ( (u0_r-u0_l) * slope < 0.0 ) then
248  u0_l = u0_r
249  u1_l = 0.0
250  u1_r = 0.0
251  endif
252 
253  ! Store edge values and slope, build cubic and check monotonicity
254  ppoly_e(i0,1) = u0_l
255  ppoly_e(i0,2) = u0_r
256  ppoly_s(i0,1) = u1_l
257  ppoly_s(i0,2) = u1_r
258 
259  ! Store edge values and slope, build cubic and check monotonicity
260  call build_cubic_interpolant( h, i0, ppoly_e, ppoly_s, ppoly_coef )
261  monotonic = is_cubic_monotonic( ppoly_coef, i0 )
262 
263  if ( .not.monotonic ) then
264  call monotonize_cubic( h0, u0_l, u0_r, 0.0, slope, slope, u1_l, u1_r )
265 
266  ! Rebuild cubic after monotonization
267  ppoly_s(i0,1) = u1_l
268  ppoly_s(i0,2) = u1_r
269  call build_cubic_interpolant( h, i0, ppoly_e, ppoly_s, ppoly_coef )
270 
271  endif
272 
273  ! ----- Right boundary -----
274  i0 = n-1
275  i1 = n
276  h0 = h(i0) + hneglect_edge
277  h1 = h(i1) + hneglect_edge
278  u0 = u(i0)
279  u1 = u(i1)
280 
281  ! Compute the right edge slope in neighboring cell and express it in
282  ! the global coordinate system
283  b = ppoly_coef(i0,2)
284  c = ppoly_coef(i0,3)
285  d = ppoly_coef(i0,4)
286  u1_l = (b + 2*c + 3*d) / ( h0 + hneglect ) ! derivative evaluated at xi = 1.0
287 
288  ! Limit the left slope by the PLM limited slope
289  slope = 2.0 * ( u1 - u0 ) / ( h1 + hneglect )
290  if ( abs(u1_l) > abs(slope) ) then
291  u1_l = slope
292  endif
293 
294  ! The left edge value in the boundary cell is taken to be the right
295  ! edge value in the neighboring cell
296  u0_l = ppoly_e(i0,2)
297 
298  ! Given the left edge value and slope, we determine the right
299  ! edge value and slope by computing the parabola as determined by
300  ! the left edge value and slope and the boundary cell average
301  u0_r = 3.0 * u1 - 0.5 * h1*u1_l - 2.0 * u0_l
302  u1_r = ( 6.0 * u1 - 2.0 * h1*u1_l - 6.0 * u0_l) / ( h1 + hneglect )
303 
304  ! Check whether the edge values are monotonic. For example, if the right edge
305  ! value is smaller than the left edge value while the slope is positive, the
306  ! edge values are inconsistent and we need to modify the right edge value
307  if ( (u0_r-u0_l) * slope < 0.0 ) then
308  u0_r = u0_l
309  u1_l = 0.0
310  u1_r = 0.0
311  endif
312 
313  ! Store edge values and slope, build cubic and check monotonicity
314  ppoly_e(i1,1) = u0_l
315  ppoly_e(i1,2) = u0_r
316  ppoly_s(i1,1) = u1_l
317  ppoly_s(i1,2) = u1_r
318 
319  call build_cubic_interpolant( h, i1, ppoly_e, ppoly_s, ppoly_coef )
320  monotonic = is_cubic_monotonic( ppoly_coef, i1 )
321 
322  if ( .not.monotonic ) then
323  call monotonize_cubic( h1, u0_l, u0_r, slope, 0.0, slope, u1_l, u1_r )
324 
325  ! Rebuild cubic after monotonization
326  ppoly_s(i1,1) = u1_l
327  ppoly_s(i1,2) = u1_r
328  call build_cubic_interpolant( h, i1, ppoly_e, ppoly_s, ppoly_coef )
329 
330  endif
331 
332 end subroutine p3m_boundary_extrapolation
333 
334 
335 !> Build cubic interpolant in cell k
336 !!
337 !! Given edge values and edge slopes, compute coefficients of cubic in cell k.
338 !!
339 !! NOTE: edge values and slopes MUST have been properly calculated prior to
340 !! calling this routine.
341 subroutine build_cubic_interpolant( h, k, ppoly_E, ppoly_S, ppoly_coef )
342  real, dimension(:), intent(in) :: h !< cell widths (size N) [H]
343  integer, intent(in) :: k !< The index of the cell to work on
344  real, dimension(:,:), intent(in) :: ppoly_E !< Edge value of polynomial in arbitrary units [A]
345  real, dimension(:,:), intent(in) :: ppoly_S !< Edge slope of polynomial [A H-1]
346  real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial [A]
347 
348  ! Local variables
349  real :: u0_l, u0_r ! edge values [A]
350  real :: u1_l, u1_r ! edge slopes times the cell width [A]
351  real :: h_c ! cell width [H]
352  real :: a0, a1, a2, a3 ! cubic coefficients [A]
353 
354  h_c = h(k)
355 
356  u0_l = ppoly_e(k,1)
357  u0_r = ppoly_e(k,2)
358 
359  u1_l = ppoly_s(k,1) * h_c
360  u1_r = ppoly_s(k,2) * h_c
361 
362  a0 = u0_l
363  a1 = u1_l
364  a2 = 3.0 * ( u0_r - u0_l ) - u1_r - 2.0 * u1_l
365  a3 = u1_r + u1_l + 2.0 * ( u0_l - u0_r )
366 
367  ppoly_coef(k,1) = a0
368  ppoly_coef(k,2) = a1
369  ppoly_coef(k,3) = a2
370  ppoly_coef(k,4) = a3
371 
372 end subroutine build_cubic_interpolant
373 
374 
375 !> Check whether the cubic reconstruction in cell k is monotonic
376 !!
377 !! This function checks whether the cubic curve in cell k is monotonic.
378 !! If so, returns 1. Otherwise, returns 0.
379 !!
380 !! The cubic is monotonic if the first derivative is single-signed in (0,1).
381 !! Hence, we check whether the roots (if any) lie inside this interval. If there
382 !! is no root or if both roots lie outside this interval, the cubic is monotonic.
383 logical function is_cubic_monotonic( ppoly_coef, k )
384  real, dimension(:,:), intent(in) :: ppoly_coef !< Coefficients of cubic polynomial in arbitary units [A]
385  integer, intent(in) :: k !< The index of the cell to work on
386  ! Local variables
387  real :: a, b, c ! Coefficients of the first derivative of the cubic [A]
388 
389  a = ppoly_coef(k,2)
390  b = 2.0 * ppoly_coef(k,3)
391  c = 3.0 * ppoly_coef(k,4)
392 
393  ! Look for real roots of the quadratic derivative equation, c*x**2 + b*x + a = 0, in (0, 1)
394  if (b*b - 4.0*a*c <= 0.0) then ! The cubic is monotonic everywhere.
395  is_cubic_monotonic = .true.
396  elseif (a * (a + (b + c)) < 0.0) then ! The derivative changes sign between the endpoints of (0, 1)
397  is_cubic_monotonic = .false.
398  elseif (b * (b + 2.0*c) < 0.0) then ! The second derivative changes sign inside of (0, 1)
399  is_cubic_monotonic = .false.
400  else
401  is_cubic_monotonic = .true.
402  endif
403 
404 end function is_cubic_monotonic
405 
406 !> Monotonize a cubic curve by modifying the edge slopes.
407 !!
408 !! This routine takes care of monotonizing a cubic on [0,1] by modifying the
409 !! edge slopes. The edge values are NOT modified. The cubic is entirely
410 !! determined by the four degrees of freedom u0_l, u0_r, u1_l and u1_r.
411 !!
412 !! u1_l and u1_r are the edge slopes expressed in the GLOBAL coordinate system.
413 !!
414 !! The monotonization occurs as follows.
415 !
416 !! 1. The edge slopes are set to 0 if they are inconsistent with the limited
417 !! PLM slope
418 !! 2. We check whether we can find an inflexion point in [0,1]. At most one
419 !! inflexion point may exist.
420 !! a. If there is no inflexion point, the cubic is monotonic.
421 !! b. If there is one inflexion point and it lies outside [0,1], the
422 !! cubic is monotonic.
423 !! c. If there is one inflexion point and it lies in [0,1] and the slope
424 !! at the location of the inflexion point is consistent, the cubic
425 !! is monotonic.
426 !! d. If the inflexion point lies in [0,1] but the slope is inconsistent,
427 !! we go to (3) to shift the location of the inflexion point to the left
428 !! or to the right. To the left when the 2nd-order left slope is smaller
429 !! than the 2nd order right slope.
430 !! 3. Edge slopes are modified to shift the inflexion point, either onto the left
431 !! edge or onto the right edge.
432 
433 subroutine monotonize_cubic( h, u0_l, u0_r, sigma_l, sigma_r, slope, u1_l, u1_r )
434  real, intent(in) :: h !< cell width [H]
435  real, intent(in) :: u0_l !< left edge value in arbitrary units [A]
436  real, intent(in) :: u0_r !< right edge value [A]
437  real, intent(in) :: sigma_l !< left 2nd-order slopes [A H-1]
438  real, intent(in) :: sigma_r !< right 2nd-order slopes [A H-1]
439  real, intent(in) :: slope !< limited PLM slope [A H-1]
440  real, intent(inout) :: u1_l !< left edge slopes [A H-1]
441  real, intent(inout) :: u1_r !< right edge slopes [A H-1]
442  ! Local variables
443  logical :: found_ip
444  logical :: inflexion_l ! bool telling if inflex. pt must be on left
445  logical :: inflexion_r ! bool telling if inflex. pt must be on right
446  real :: a1, a2, a3 ! Temporary slopes times the cell width [A]
447  real :: u1_l_tmp ! trial left edge slope [A H-1]
448  real :: u1_r_tmp ! trial right edge slope [A H-1]
449  real :: xi_ip ! location of inflexion point in cell coordinates (0,1) [nondim]
450  real :: slope_ip ! slope at inflexion point times cell width [A]
451 
452  found_ip = .false.
453  inflexion_l = .false.
454  inflexion_r = .false.
455 
456  ! If the edge slopes are inconsistent w.r.t. the limited PLM slope,
457  ! set them to zero
458  if ( u1_l*slope <= 0.0 ) then
459  u1_l = 0.0
460  endif
461 
462  if ( u1_r*slope <= 0.0 ) then
463  u1_r = 0.0
464  endif
465 
466  ! Compute the location of the inflexion point, which is the root
467  ! of the second derivative
468  a1 = h * u1_l
469  a2 = 3.0 * ( u0_r - u0_l ) - h*(u1_r + 2.0*u1_l)
470  a3 = h*(u1_r + u1_l) + 2.0*(u0_l - u0_r)
471 
472  ! There is a possible root (and inflexion point) only if a3 is nonzero.
473  ! When a3 is zero, the second derivative of the cubic is constant (the
474  ! cubic degenerates into a parabola) and no inflexion point exists.
475  if ( a3 /= 0.0 ) then
476  ! Location of inflexion point
477  xi_ip = - a2 / (3.0 * a3)
478 
479  ! If the inflexion point lies in [0,1], change boolean value
480  if ( (xi_ip >= 0.0) .AND. (xi_ip <= 1.0) ) then
481  found_ip = .true.
482  endif
483  endif
484 
485  ! When there is an inflexion point within [0,1], check the slope
486  ! to see if it is consistent with the limited PLM slope. If not,
487  ! decide on which side we want to collapse the inflexion point.
488  ! If the inflexion point lies on one of the edges, the cubic is
489  ! guaranteed to be monotonic
490  if ( found_ip ) then
491  slope_ip = a1 + 2.0*a2*xi_ip + 3.0*a3*xi_ip*xi_ip
492 
493  ! Check whether slope is consistent
494  if ( slope_ip*slope < 0.0 ) then
495  if ( abs(sigma_l) < abs(sigma_r) ) then
496  inflexion_l = .true.
497  else
498  inflexion_r = .true.
499  endif
500  endif
501  endif ! found_ip
502 
503  ! At this point, if the cubic is not monotonic, we know where the
504  ! inflexion point should lie. When the cubic is monotonic, both
505  ! 'inflexion_l' and 'inflexion_r' are false and nothing is to be done.
506 
507  ! Move inflexion point on the left
508  if ( inflexion_l ) then
509 
510  u1_l_tmp = 1.5*(u0_r-u0_l)/h - 0.5*u1_r
511  u1_r_tmp = 3.0*(u0_r-u0_l)/h - 2.0*u1_l
512 
513  if ( (u1_l_tmp*slope < 0.0) .AND. (u1_r_tmp*slope < 0.0) ) then
514 
515  u1_l = 0.0
516  u1_r = 3.0 * (u0_r - u0_l) / h
517 
518  elseif (u1_l_tmp*slope < 0.0) then
519 
520  u1_r = u1_r_tmp
521  u1_l = 1.5*(u0_r - u0_l)/h - 0.5*u1_r
522 
523  elseif (u1_r_tmp*slope < 0.0) then
524 
525  u1_l = u1_l_tmp
526  u1_r = 3.0*(u0_r - u0_l)/h - 2.0*u1_l
527 
528  else
529 
530  u1_l = u1_l_tmp
531  u1_r = u1_r_tmp
532 
533  endif
534 
535  endif ! end treating case with inflexion point on the left
536 
537  ! Move inflexion point on the right
538  if ( inflexion_r ) then
539 
540  u1_l_tmp = 3.0*(u0_r-u0_l)/h - 2.0*u1_r
541  u1_r_tmp = 1.5*(u0_r-u0_l)/h - 0.5*u1_l
542 
543  if ( (u1_l_tmp*slope < 0.0) .AND. (u1_r_tmp*slope < 0.0) ) then
544 
545  u1_l = 3.0 * (u0_r - u0_l) / h
546  u1_r = 0.0
547 
548  elseif (u1_l_tmp*slope < 0.0) then
549 
550  u1_r = u1_r_tmp
551  u1_l = 3.0*(u0_r - u0_l)/h - 2.0*u1_r
552 
553  elseif (u1_r_tmp*slope < 0.0) then
554 
555  u1_l = u1_l_tmp
556  u1_r = 1.5*(u0_r - u0_l)/h - 0.5*u1_l
557 
558  else
559 
560  u1_l = u1_l_tmp
561  u1_r = u1_r_tmp
562 
563  endif
564 
565  endif ! end treating case with inflexion point on the right
566 
567  ! Zero out negligibly small slopes.
568  if ( abs(u1_l*h) < epsilon(u0_l) * (abs(u0_l) + abs(u0_r)) ) u1_l = 0.0
569  if ( abs(u1_r*h) < epsilon(u0_l) * (abs(u0_l) + abs(u0_r)) ) u1_r = 0.0
570 
571 end subroutine monotonize_cubic
572 
573 !> \namespace p3m_functions
574 !!
575 !! Date of creation: 2008.06.09
576 !! L. White
577 !!
578 !! This module contains p3m interpolation routines.
579 !!
580 !! p3m interpolation is performed by estimating the edge values and slopes
581 !! and constructing a cubic polynomial. We then make sure that the edge values
582 !! are bounded and continuous and we then modify the slopes to get a monotonic
583 !! cubic curve.
584 
585 end module p3m_functions
p3m_functions::build_cubic_interpolant
subroutine build_cubic_interpolant(h, k, ppoly_E, ppoly_S, ppoly_coef)
Build cubic interpolant in cell k.
Definition: P3M_functions.F90:342
p3m_functions::is_cubic_monotonic
logical function is_cubic_monotonic(ppoly_coef, k)
Check whether the cubic reconstruction in cell k is monotonic.
Definition: P3M_functions.F90:384
p3m_functions::p3m_interpolation
subroutine, public p3m_interpolation(N, h, u, ppoly_E, ppoly_S, ppoly_coef, h_neglect)
Set up a piecewise cubic interpolation from cell averages and estimated edge slopes and values.
Definition: P3M_functions.F90:29
p3m_functions
Cubic interpolation functions.
Definition: P3M_functions.F90:2
p3m_functions::monotonize_cubic
subroutine monotonize_cubic(h, u0_l, u0_r, sigma_l, sigma_r, slope, u1_l, u1_r)
Monotonize a cubic curve by modifying the edge slopes.
Definition: P3M_functions.F90:434
regrid_edge_values::average_discontinuous_edge_values
subroutine, public average_discontinuous_edge_values(N, edge_val)
Replace discontinuous collocated edge values with their average.
Definition: regrid_edge_values.F90:141
regrid_edge_values
Edge value estimation for high-order resconstruction.
Definition: regrid_edge_values.F90:2
p3m_functions::hneglect_edge_dflt
real, parameter hneglect_edge_dflt
Default value of a negligible edge thickness.
Definition: P3M_functions.F90:14
p3m_functions::p3m_limiter
subroutine p3m_limiter(N, h, u, ppoly_E, ppoly_S, ppoly_coef, h_neglect)
Adust a piecewise cubic reconstruction with a limiter that adjusts the edge values and slopes.
Definition: P3M_functions.F90:61
regrid_edge_values::bound_edge_values
subroutine, public bound_edge_values(N, h, u, edge_val, h_neglect)
Bound edge values by neighboring cell averages.
Definition: regrid_edge_values.F90:48
p3m_functions::hneglect_dflt
real, parameter hneglect_dflt
Default value of a negligible cell thickness.
Definition: P3M_functions.F90:13
p3m_functions::p3m_boundary_extrapolation
subroutine, public p3m_boundary_extrapolation(N, h, u, ppoly_E, ppoly_S, ppoly_coef, h_neglect, h_neglect_edge)
Calculate the edge values and slopes at boundary cells as part of building a piecewise cubic sub-grid...
Definition: P3M_functions.F90:190