MOM6
Functions/Subroutines | Variables
p3m_functions Module Reference

Detailed Description

Cubic interpolation functions.

Date of creation: 2008.06.09 L. White.

This module contains p3m interpolation routines.

p3m interpolation is performed by estimating the edge values and slopes and constructing a cubic polynomial. We then make sure that the edge values are bounded and continuous and we then modify the slopes to get a monotonic cubic curve.

Functions/Subroutines

subroutine, public p3m_interpolation (N, h, u, edge_values, ppoly_S, ppoly_coef, h_neglect, answers_2018)
 Set up a piecewise cubic interpolation from cell averages and estimated edge slopes and values. More...
 
subroutine p3m_limiter (N, h, u, edge_values, ppoly_S, ppoly_coef, h_neglect, answers_2018)
 Adust a piecewise cubic reconstruction with a limiter that adjusts the edge values and slopes. More...
 
subroutine, public p3m_boundary_extrapolation (N, h, u, edge_values, 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 scale profiles. More...
 
subroutine build_cubic_interpolant (h, k, edge_values, ppoly_S, ppoly_coef)
 Build cubic interpolant in cell k. More...
 
logical function is_cubic_monotonic (ppoly_coef, k)
 Check whether the cubic reconstruction in cell k is monotonic. More...
 
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. More...
 

Variables

real, parameter hneglect_dflt = 1.E-30
 Default value of a negligible cell thickness.
 
real, parameter hneglect_edge_dflt = 1.E-10
 Default value of a negligible edge thickness.
 

Function/Subroutine Documentation

◆ build_cubic_interpolant()

subroutine p3m_functions::build_cubic_interpolant ( real, dimension(:), intent(in)  h,
integer, intent(in)  k,
real, dimension(:,:), intent(in)  edge_values,
real, dimension(:,:), intent(in)  ppoly_S,
real, dimension(:,:), intent(inout)  ppoly_coef 
)
private

Build cubic interpolant in cell k.

Given edge values and edge slopes, compute coefficients of cubic in cell k.

NOTE: edge values and slopes MUST have been properly calculated prior to calling this routine.

Parameters
[in]hcell widths (size N) [H]
[in]kThe index of the cell to work on
[in]edge_valuesEdge value of polynomial in arbitrary units [A]
[in]ppoly_sEdge slope of polynomial [A H-1]
[in,out]ppoly_coefCoefficients of polynomial [A]

Definition at line 344 of file P3M_functions.F90.

344  real, dimension(:), intent(in) :: h !< cell widths (size N) [H]
345  integer, intent(in) :: k !< The index of the cell to work on
346  real, dimension(:,:), intent(in) :: edge_values !< Edge value of polynomial in arbitrary units [A]
347  real, dimension(:,:), intent(in) :: ppoly_S !< Edge slope of polynomial [A H-1]
348  real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial [A]
349 
350  ! Local variables
351  real :: u0_l, u0_r ! edge values [A]
352  real :: u1_l, u1_r ! edge slopes times the cell width [A]
353  real :: h_c ! cell width [H]
354  real :: a0, a1, a2, a3 ! cubic coefficients [A]
355 
356  h_c = h(k)
357 
358  u0_l = edge_values(k,1)
359  u0_r = edge_values(k,2)
360 
361  u1_l = ppoly_s(k,1) * h_c
362  u1_r = ppoly_s(k,2) * h_c
363 
364  a0 = u0_l
365  a1 = u1_l
366  a2 = 3.0 * ( u0_r - u0_l ) - u1_r - 2.0 * u1_l
367  a3 = u1_r + u1_l + 2.0 * ( u0_l - u0_r )
368 
369  ppoly_coef(k,1) = a0
370  ppoly_coef(k,2) = a1
371  ppoly_coef(k,3) = a2
372  ppoly_coef(k,4) = a3
373 

◆ is_cubic_monotonic()

logical function p3m_functions::is_cubic_monotonic ( real, dimension(:,:), intent(in)  ppoly_coef,
integer, intent(in)  k 
)
private

Check whether the cubic reconstruction in cell k is monotonic.

This function checks whether the cubic curve in cell k is monotonic. If so, returns 1. Otherwise, returns 0.

The cubic is monotonic if the first derivative is single-signed in (0,1). Hence, we check whether the roots (if any) lie inside this interval. If there is no root or if both roots lie outside this interval, the cubic is monotonic.

Parameters
[in]ppoly_coefCoefficients of cubic polynomial in arbitary units [A]
[in]kThe index of the cell to work on

Definition at line 386 of file P3M_functions.F90.

386  real, dimension(:,:), intent(in) :: ppoly_coef !< Coefficients of cubic polynomial in arbitary units [A]
387  integer, intent(in) :: k !< The index of the cell to work on
388  ! Local variables
389  real :: a, b, c ! Coefficients of the first derivative of the cubic [A]
390 
391  a = ppoly_coef(k,2)
392  b = 2.0 * ppoly_coef(k,3)
393  c = 3.0 * ppoly_coef(k,4)
394 
395  ! Look for real roots of the quadratic derivative equation, c*x**2 + b*x + a = 0, in (0, 1)
396  if (b*b - 4.0*a*c <= 0.0) then ! The cubic is monotonic everywhere.
397  is_cubic_monotonic = .true.
398  elseif (a * (a + (b + c)) < 0.0) then ! The derivative changes sign between the endpoints of (0, 1)
399  is_cubic_monotonic = .false.
400  elseif (b * (b + 2.0*c) < 0.0) then ! The second derivative changes sign inside of (0, 1)
401  is_cubic_monotonic = .false.
402  else
403  is_cubic_monotonic = .true.
404  endif
405 

◆ monotonize_cubic()

subroutine p3m_functions::monotonize_cubic ( real, intent(in)  h,
real, intent(in)  u0_l,
real, intent(in)  u0_r,
real, intent(in)  sigma_l,
real, intent(in)  sigma_r,
real, intent(in)  slope,
real, intent(inout)  u1_l,
real, intent(inout)  u1_r 
)
private

Monotonize a cubic curve by modifying the edge slopes.

This routine takes care of monotonizing a cubic on [0,1] by modifying the edge slopes. The edge values are NOT modified. The cubic is entirely determined by the four degrees of freedom u0_l, u0_r, u1_l and u1_r.

u1_l and u1_r are the edge slopes expressed in the GLOBAL coordinate system.

The monotonization occurs as follows.

Parameters
[in]hcell width [H]
[in]u0_lleft edge value in arbitrary units [A]
[in]u0_rright edge value [A]
[in]sigma_lleft 2nd-order slopes [A H-1]
[in]sigma_rright 2nd-order slopes [A H-1]
[in]slopelimited PLM slope [A H-1]
[in,out]u1_lleft edge slopes [A H-1]
[in,out]u1_rright edge slopes [A H-1]

Definition at line 436 of file P3M_functions.F90.

436  real, intent(in) :: h !< cell width [H]
437  real, intent(in) :: u0_l !< left edge value in arbitrary units [A]
438  real, intent(in) :: u0_r !< right edge value [A]
439  real, intent(in) :: sigma_l !< left 2nd-order slopes [A H-1]
440  real, intent(in) :: sigma_r !< right 2nd-order slopes [A H-1]
441  real, intent(in) :: slope !< limited PLM slope [A H-1]
442  real, intent(inout) :: u1_l !< left edge slopes [A H-1]
443  real, intent(inout) :: u1_r !< right edge slopes [A H-1]
444  ! Local variables
445  logical :: found_ip
446  logical :: inflexion_l ! bool telling if inflex. pt must be on left
447  logical :: inflexion_r ! bool telling if inflex. pt must be on right
448  real :: a1, a2, a3 ! Temporary slopes times the cell width [A]
449  real :: u1_l_tmp ! trial left edge slope [A H-1]
450  real :: u1_r_tmp ! trial right edge slope [A H-1]
451  real :: xi_ip ! location of inflexion point in cell coordinates (0,1) [nondim]
452  real :: slope_ip ! slope at inflexion point times cell width [A]
453 
454  found_ip = .false.
455  inflexion_l = .false.
456  inflexion_r = .false.
457 
458  ! If the edge slopes are inconsistent w.r.t. the limited PLM slope,
459  ! set them to zero
460  if ( u1_l*slope <= 0.0 ) then
461  u1_l = 0.0
462  endif
463 
464  if ( u1_r*slope <= 0.0 ) then
465  u1_r = 0.0
466  endif
467 
468  ! Compute the location of the inflexion point, which is the root
469  ! of the second derivative
470  a1 = h * u1_l
471  a2 = 3.0 * ( u0_r - u0_l ) - h*(u1_r + 2.0*u1_l)
472  a3 = h*(u1_r + u1_l) + 2.0*(u0_l - u0_r)
473 
474  ! There is a possible root (and inflexion point) only if a3 is nonzero.
475  ! When a3 is zero, the second derivative of the cubic is constant (the
476  ! cubic degenerates into a parabola) and no inflexion point exists.
477  if ( a3 /= 0.0 ) then
478  ! Location of inflexion point
479  xi_ip = - a2 / (3.0 * a3)
480 
481  ! If the inflexion point lies in [0,1], change boolean value
482  if ( (xi_ip >= 0.0) .AND. (xi_ip <= 1.0) ) then
483  found_ip = .true.
484  endif
485  endif
486 
487  ! When there is an inflexion point within [0,1], check the slope
488  ! to see if it is consistent with the limited PLM slope. If not,
489  ! decide on which side we want to collapse the inflexion point.
490  ! If the inflexion point lies on one of the edges, the cubic is
491  ! guaranteed to be monotonic
492  if ( found_ip ) then
493  slope_ip = a1 + 2.0*a2*xi_ip + 3.0*a3*xi_ip*xi_ip
494 
495  ! Check whether slope is consistent
496  if ( slope_ip*slope < 0.0 ) then
497  if ( abs(sigma_l) < abs(sigma_r) ) then
498  inflexion_l = .true.
499  else
500  inflexion_r = .true.
501  endif
502  endif
503  endif ! found_ip
504 
505  ! At this point, if the cubic is not monotonic, we know where the
506  ! inflexion point should lie. When the cubic is monotonic, both
507  ! 'inflexion_l' and 'inflexion_r' are false and nothing is to be done.
508 
509  ! Move inflexion point on the left
510  if ( inflexion_l ) then
511 
512  u1_l_tmp = 1.5*(u0_r-u0_l)/h - 0.5*u1_r
513  u1_r_tmp = 3.0*(u0_r-u0_l)/h - 2.0*u1_l
514 
515  if ( (u1_l_tmp*slope < 0.0) .AND. (u1_r_tmp*slope < 0.0) ) then
516 
517  u1_l = 0.0
518  u1_r = 3.0 * (u0_r - u0_l) / h
519 
520  elseif (u1_l_tmp*slope < 0.0) then
521 
522  u1_r = u1_r_tmp
523  u1_l = 1.5*(u0_r - u0_l)/h - 0.5*u1_r
524 
525  elseif (u1_r_tmp*slope < 0.0) then
526 
527  u1_l = u1_l_tmp
528  u1_r = 3.0*(u0_r - u0_l)/h - 2.0*u1_l
529 
530  else
531 
532  u1_l = u1_l_tmp
533  u1_r = u1_r_tmp
534 
535  endif
536 
537  endif ! end treating case with inflexion point on the left
538 
539  ! Move inflexion point on the right
540  if ( inflexion_r ) then
541 
542  u1_l_tmp = 3.0*(u0_r-u0_l)/h - 2.0*u1_r
543  u1_r_tmp = 1.5*(u0_r-u0_l)/h - 0.5*u1_l
544 
545  if ( (u1_l_tmp*slope < 0.0) .AND. (u1_r_tmp*slope < 0.0) ) then
546 
547  u1_l = 3.0 * (u0_r - u0_l) / h
548  u1_r = 0.0
549 
550  elseif (u1_l_tmp*slope < 0.0) then
551 
552  u1_r = u1_r_tmp
553  u1_l = 3.0*(u0_r - u0_l)/h - 2.0*u1_r
554 
555  elseif (u1_r_tmp*slope < 0.0) then
556 
557  u1_l = u1_l_tmp
558  u1_r = 1.5*(u0_r - u0_l)/h - 0.5*u1_l
559 
560  else
561 
562  u1_l = u1_l_tmp
563  u1_r = u1_r_tmp
564 
565  endif
566 
567  endif ! end treating case with inflexion point on the right
568 
569  ! Zero out negligibly small slopes.
570  if ( abs(u1_l*h) < epsilon(u0_l) * (abs(u0_l) + abs(u0_r)) ) u1_l = 0.0
571  if ( abs(u1_r*h) < epsilon(u0_l) * (abs(u0_l) + abs(u0_r)) ) u1_r = 0.0
572 

◆ p3m_boundary_extrapolation()

subroutine, public p3m_functions::p3m_boundary_extrapolation ( integer, intent(in)  N,
real, dimension(:), intent(in)  h,
real, dimension(:), intent(in)  u,
real, dimension(:,:), intent(inout)  edge_values,
real, dimension(:,:), intent(inout)  ppoly_S,
real, dimension(:,:), intent(inout)  ppoly_coef,
real, intent(in), optional  h_neglect,
real, intent(in), optional  h_neglect_edge 
)

Calculate the edge values and slopes at boundary cells as part of building a piecewise cubic sub-grid scale profiles.

The following explanations apply to the left boundary cell. The same reasoning holds for the right boundary cell.

A cubic needs to be built in the cell and requires four degrees of freedom, which are the edge values and slopes. The right edge values and slopes are taken to be that of the neighboring cell (i.e., the left edge value and slope of the neighboring cell). The left edge value and slope are determined by computing the parabola based on the cell average and the right edge value and slope. The resulting cubic is not necessarily monotonic and the slopes are subsequently modified to yield a monotonic cubic.

Parameters
[in]nNumber of cells
[in]hcell widths (size N) [H]
[in]ucell averages (size N) in arbitrary units [A]
[in,out]edge_valuesEdge value of polynomial [A]
[in,out]ppoly_sEdge slope of polynomial [A H-1]
[in,out]ppoly_coefCoefficients of polynomial [A]
[in]h_neglectA negligibly small width for the purpose of cell reconstructions [H]
[in]h_neglect_edgeA negligibly small width for the purpose of finding edge values [H]

Definition at line 193 of file P3M_functions.F90.

193  integer, intent(in) :: N !< Number of cells
194  real, dimension(:), intent(in) :: h !< cell widths (size N) [H]
195  real, dimension(:), intent(in) :: u !< cell averages (size N) in arbitrary units [A]
196  real, dimension(:,:), intent(inout) :: edge_values !< Edge value of polynomial [A]
197  real, dimension(:,:), intent(inout) :: ppoly_S !< Edge slope of polynomial [A H-1]
198  real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial [A]
199  real, optional, intent(in) :: h_neglect !< A negligibly small width for the
200  !! purpose of cell reconstructions [H]
201  real, optional, intent(in) :: h_neglect_edge !< A negligibly small width
202  !! for the purpose of finding edge values [H]
203  ! Local variables
204  integer :: i0, i1
205  logical :: monotonic ! boolean indicating whether the cubic is monotonic
206  real :: u0, u1 ! Values of u in two adjacent cells [A]
207  real :: h0, h1 ! Values of h in two adjacent cells, plus a smal increment [H]
208  real :: b, c, d ! Temporary variables [A]
209  real :: u0_l, u0_r ! Left and right edge values [A]
210  real :: u1_l, u1_r ! Left and right edge slopes [A H-1]
211  real :: slope ! The cell center slope [A H-1]
212  real :: hNeglect, hNeglect_edge ! Negligibly small thickness [H]
213 
214  hneglect = hneglect_dflt ; if (present(h_neglect)) hneglect = h_neglect
215  hneglect_edge = hneglect_edge_dflt ; if (present(h_neglect_edge)) hneglect_edge = h_neglect_edge
216 
217  ! ----- Left boundary -----
218  i0 = 1
219  i1 = 2
220  h0 = h(i0) + hneglect_edge
221  h1 = h(i1) + hneglect_edge
222  u0 = u(i0)
223  u1 = u(i1)
224 
225  ! Compute the left edge slope in neighboring cell and express it in
226  ! the global coordinate system
227  b = ppoly_coef(i1,2)
228  u1_r = b / h1 ! derivative evaluated at xi = 0.0, expressed w.r.t. x
229 
230  ! Limit the right slope by the PLM limited slope
231  slope = 2.0 * ( u1 - u0 ) / ( h0 + hneglect )
232  if ( abs(u1_r) > abs(slope) ) then
233  u1_r = slope
234  endif
235 
236  ! The right edge value in the boundary cell is taken to be the left
237  ! edge value in the neighboring cell
238  u0_r = edge_values(i1,1)
239 
240  ! Given the right edge value and slope, we determine the left
241  ! edge value and slope by computing the parabola as determined by
242  ! the right edge value and slope and the boundary cell average
243  u0_l = 3.0 * u0 + 0.5 * h0*u1_r - 2.0 * u0_r
244  u1_l = ( - 6.0 * u0 - 2.0 * h0*u1_r + 6.0 * u0_r) / ( h0 + hneglect )
245 
246  ! Check whether the edge values are monotonic. For example, if the left edge
247  ! value is larger than the right edge value while the slope is positive, the
248  ! edge values are inconsistent and we need to modify the left edge value
249  if ( (u0_r-u0_l) * slope < 0.0 ) then
250  u0_l = u0_r
251  u1_l = 0.0
252  u1_r = 0.0
253  endif
254 
255  ! Store edge values and slope, build cubic and check monotonicity
256  edge_values(i0,1) = u0_l
257  edge_values(i0,2) = u0_r
258  ppoly_s(i0,1) = u1_l
259  ppoly_s(i0,2) = u1_r
260 
261  ! Store edge values and slope, build cubic and check monotonicity
262  call build_cubic_interpolant( h, i0, edge_values, ppoly_s, ppoly_coef )
263  monotonic = is_cubic_monotonic( ppoly_coef, i0 )
264 
265  if ( .not.monotonic ) then
266  call monotonize_cubic( h0, u0_l, u0_r, 0.0, slope, slope, u1_l, u1_r )
267 
268  ! Rebuild cubic after monotonization
269  ppoly_s(i0,1) = u1_l
270  ppoly_s(i0,2) = u1_r
271  call build_cubic_interpolant( h, i0, edge_values, ppoly_s, ppoly_coef )
272 
273  endif
274 
275  ! ----- Right boundary -----
276  i0 = n-1
277  i1 = n
278  h0 = h(i0) + hneglect_edge
279  h1 = h(i1) + hneglect_edge
280  u0 = u(i0)
281  u1 = u(i1)
282 
283  ! Compute the right edge slope in neighboring cell and express it in
284  ! the global coordinate system
285  b = ppoly_coef(i0,2)
286  c = ppoly_coef(i0,3)
287  d = ppoly_coef(i0,4)
288  u1_l = (b + 2*c + 3*d) / ( h0 + hneglect ) ! derivative evaluated at xi = 1.0
289 
290  ! Limit the left slope by the PLM limited slope
291  slope = 2.0 * ( u1 - u0 ) / ( h1 + hneglect )
292  if ( abs(u1_l) > abs(slope) ) then
293  u1_l = slope
294  endif
295 
296  ! The left edge value in the boundary cell is taken to be the right
297  ! edge value in the neighboring cell
298  u0_l = edge_values(i0,2)
299 
300  ! Given the left edge value and slope, we determine the right
301  ! edge value and slope by computing the parabola as determined by
302  ! the left edge value and slope and the boundary cell average
303  u0_r = 3.0 * u1 - 0.5 * h1*u1_l - 2.0 * u0_l
304  u1_r = ( 6.0 * u1 - 2.0 * h1*u1_l - 6.0 * u0_l) / ( h1 + hneglect )
305 
306  ! Check whether the edge values are monotonic. For example, if the right edge
307  ! value is smaller than the left edge value while the slope is positive, the
308  ! edge values are inconsistent and we need to modify the right edge value
309  if ( (u0_r-u0_l) * slope < 0.0 ) then
310  u0_r = u0_l
311  u1_l = 0.0
312  u1_r = 0.0
313  endif
314 
315  ! Store edge values and slope, build cubic and check monotonicity
316  edge_values(i1,1) = u0_l
317  edge_values(i1,2) = u0_r
318  ppoly_s(i1,1) = u1_l
319  ppoly_s(i1,2) = u1_r
320 
321  call build_cubic_interpolant( h, i1, edge_values, ppoly_s, ppoly_coef )
322  monotonic = is_cubic_monotonic( ppoly_coef, i1 )
323 
324  if ( .not.monotonic ) then
325  call monotonize_cubic( h1, u0_l, u0_r, slope, 0.0, slope, u1_l, u1_r )
326 
327  ! Rebuild cubic after monotonization
328  ppoly_s(i1,1) = u1_l
329  ppoly_s(i1,2) = u1_r
330  call build_cubic_interpolant( h, i1, edge_values, ppoly_s, ppoly_coef )
331 
332  endif
333 

◆ p3m_interpolation()

subroutine, public p3m_functions::p3m_interpolation ( integer, intent(in)  N,
real, dimension(:), intent(in)  h,
real, dimension(:), intent(in)  u,
real, dimension(:,:), intent(inout)  edge_values,
real, dimension(:,:), intent(inout)  ppoly_S,
real, dimension(:,:), intent(inout)  ppoly_coef,
real, intent(in), optional  h_neglect,
logical, intent(in), optional  answers_2018 
)

Set up a piecewise cubic interpolation from cell averages and estimated edge slopes and values.

Cubic interpolation between edges.

The edge values and slopes MUST have been estimated prior to calling this routine.

It is assumed that the size of the array 'u' is equal to the number of cells defining 'grid' and 'ppoly'. No consistency check is performed here.

Parameters
[in]nNumber of cells
[in]hcell widths (size N) [H]
[in]ucell averages (size N) in arbitrary units [A]
[in,out]edge_valuesEdge value of polynomial [A]
[in,out]ppoly_sEdge slope of polynomial [A H-1].
[in,out]ppoly_coefCoefficients of polynomial [A]
[in]h_neglectA negligibly small width for the purpose of cell reconstructions [H]
[in]answers_2018If true use older, less acccurate expressions.

Definition at line 29 of file P3M_functions.F90.

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) :: edge_values !< 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  logical, optional, intent(in) :: answers_2018 !< If true use older, less acccurate expressions.
38 
39  ! Call the limiter for p3m, which takes care of everything from
40  ! computing the coefficients of the cubic to monotonizing it.
41  ! This routine could be called directly instead of having to call
42  ! 'P3M_interpolation' first but we do that to provide an homogeneous
43  ! interface.
44  call p3m_limiter( n, h, u, edge_values, ppoly_s, ppoly_coef, h_neglect, answers_2018 )
45 

◆ p3m_limiter()

subroutine p3m_functions::p3m_limiter ( integer, intent(in)  N,
real, dimension(:), intent(in)  h,
real, dimension(:), intent(in)  u,
real, dimension(:,:), intent(inout)  edge_values,
real, dimension(:,:), intent(inout)  ppoly_S,
real, dimension(:,:), intent(inout)  ppoly_coef,
real, intent(in), optional  h_neglect,
logical, intent(in), optional  answers_2018 
)
private

Adust a piecewise cubic reconstruction with a limiter that adjusts the edge values and slopes.

The p3m limiter operates as follows:

  1. Edge values are bounded
  2. Discontinuous edge values are systematically averaged
  3. Loop on cells and do the following a. Build cubic curve b. Check if cubic curve is monotonic c. If not, monotonize cubic curve and rebuild it

Step 3 of the monotonization process leaves all edge values unchanged.

Parameters
[in]nNumber of cells
[in]hcell widths (size N) [H]
[in]ucell averages (size N) in arbitrary units [A]
[in,out]edge_valuesEdge value of polynomial [A]
[in,out]ppoly_sEdge slope of polynomial [A H-1]
[in,out]ppoly_coefCoefficients of polynomial [A]
[in]h_neglectA negligibly small width for the purpose of cell reconstructions [H]
[in]answers_2018If true use older, less acccurate expressions.

Definition at line 62 of file P3M_functions.F90.

62  integer, intent(in) :: N !< Number of cells
63  real, dimension(:), intent(in) :: h !< cell widths (size N) [H]
64  real, dimension(:), intent(in) :: u !< cell averages (size N) in arbitrary units [A]
65  real, dimension(:,:), intent(inout) :: edge_values !< Edge value of polynomial [A]
66  real, dimension(:,:), intent(inout) :: ppoly_S !< Edge slope of polynomial [A H-1]
67  real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial [A]
68  real, optional, intent(in) :: h_neglect !< A negligibly small width for
69  !! the purpose of cell reconstructions [H]
70  logical, optional, intent(in) :: answers_2018 !< If true use older, less acccurate expressions.
71 
72  ! Local variables
73  integer :: k ! loop index
74  logical :: monotonic ! boolean indicating whether the cubic is monotonic
75  real :: u0_l, u0_r ! edge values [A]
76  real :: u1_l, u1_r ! edge slopes [A H-1]
77  real :: u_l, u_c, u_r ! left, center and right cell averages [A]
78  real :: h_l, h_c, h_r ! left, center and right cell widths [H]
79  real :: sigma_l, sigma_c, sigma_r ! left, center and right van Leer slopes [A H-1]
80  real :: slope ! retained PLM slope [A H-1]
81  real :: eps
82  real :: hNeglect ! A negligibly small thickness [H]
83 
84  hneglect = hneglect_dflt ; if (present(h_neglect)) hneglect = h_neglect
85 
86  eps = 1e-10
87 
88  ! 1. Bound edge values (boundary cells are assumed to be local extrema)
89  call bound_edge_values( n, h, u, edge_values, hneglect, answers_2018 )
90 
91  ! 2. Systematically average discontinuous edge values
92  call average_discontinuous_edge_values( n, edge_values )
93 
94 
95  ! 3. Loop on cells and do the following
96  ! (a) Build cubic curve
97  ! (b) Check if cubic curve is monotonic
98  ! (c) If not, monotonize cubic curve and rebuild it
99  do k = 1,n
100 
101  ! Get edge values, edge slopes and cell width
102  u0_l = edge_values(k,1)
103  u0_r = edge_values(k,2)
104  u1_l = ppoly_s(k,1)
105  u1_r = ppoly_s(k,2)
106 
107  ! Get cell widths and cell averages (boundary cells are assumed to
108  ! be local extrema for the sake of slopes)
109  u_c = u(k)
110  h_c = h(k)
111 
112  if ( k == 1 ) then
113  h_l = h(k)
114  u_l = u(k)
115  else
116  h_l = h(k-1)
117  u_l = u(k-1)
118  endif
119 
120  if ( k == n ) then
121  h_r = h(k)
122  u_r = u(k)
123  else
124  h_r = h(k+1)
125  u_r = u(k+1)
126  endif
127 
128  ! Compute limited slope
129  sigma_l = 2.0 * ( u_c - u_l ) / ( h_c + hneglect )
130  sigma_c = 2.0 * ( u_r - u_l ) / ( h_l + 2.0*h_c + h_r + hneglect )
131  sigma_r = 2.0 * ( u_r - u_c ) / ( h_c + hneglect )
132 
133  if ( (sigma_l * sigma_r) > 0.0 ) then
134  slope = sign( min(abs(sigma_l),abs(sigma_c),abs(sigma_r)), sigma_c )
135  else
136  slope = 0.0
137  endif
138 
139  ! If the slopes are small, set them to zero to prevent asymmetric representation near extrema.
140  if ( abs(u1_l*h_c) < epsilon(u_c)*abs(u_c) ) u1_l = 0.0
141  if ( abs(u1_r*h_c) < epsilon(u_c)*abs(u_c) ) u1_r = 0.0
142 
143  ! The edge slopes are limited from above by the respective
144  ! one-sided slopes
145  if ( abs(u1_l) > abs(sigma_l) ) then
146  u1_l = sigma_l
147  endif
148 
149  if ( abs(u1_r) > abs(sigma_r) ) then
150  u1_r = sigma_r
151  endif
152 
153  ! Build cubic interpolant (compute the coefficients)
154  call build_cubic_interpolant( h, k, edge_values, ppoly_s, ppoly_coef )
155 
156  ! Check whether cubic is monotonic
157  monotonic = is_cubic_monotonic( ppoly_coef, k )
158 
159  ! If cubic is not monotonic, monotonize it by modifiying the
160  ! edge slopes, store the new edge slopes and recompute the
161  ! cubic coefficients
162  if ( .not.monotonic ) then
163  call monotonize_cubic( h_c, u0_l, u0_r, sigma_l, sigma_r, slope, u1_l, u1_r )
164  endif
165 
166  ! Store edge slopes
167  ppoly_s(k,1) = u1_l
168  ppoly_s(k,2) = u1_r
169 
170  ! Recompute coefficients of cubic
171  call build_cubic_interpolant( h, k, edge_values, ppoly_s, ppoly_coef )
172 
173  enddo ! loop on cells
174