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, ppoly_E, ppoly_S, ppoly_coef, h_neglect)
 Set up a piecewise cubic interpolation from cell averages and estimated edge slopes and values. More...
 
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. More...
 
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 scale profiles. More...
 
subroutine build_cubic_interpolant (h, k, ppoly_E, 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. More...
 
real, parameter hneglect_edge_dflt = 1.E-10
 Default value of a negligible edge thickness. More...
 

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)  ppoly_E,
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]ppoly_eEdge 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 342 of file P3M_functions.F90.

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 

Referenced by p3m_boundary_extrapolation(), and p3m_limiter().

Here is the caller graph for this function:

◆ 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 384 of file P3M_functions.F90.

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 

Referenced by p3m_boundary_extrapolation(), and p3m_limiter().

Here is the caller graph for this function:

◆ 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 434 of file P3M_functions.F90.

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 

Referenced by p3m_boundary_extrapolation(), and p3m_limiter().

Here is the caller graph for this function:

◆ 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)  ppoly_E,
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]ppoly_eEdge 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 190 of file P3M_functions.F90.

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 

References build_cubic_interpolant(), hneglect_dflt, hneglect_edge_dflt, is_cubic_monotonic(), and monotonize_cubic().

Referenced by regrid_interp::regridding_set_ppolys().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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)  ppoly_E,
real, dimension(:,:), intent(inout)  ppoly_S,
real, dimension(:,:), intent(inout)  ppoly_coef,
real, intent(in), optional  h_neglect 
)

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]ppoly_eEdge 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]

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) :: 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 

References p3m_limiter().

Referenced by regrid_interp::regridding_set_ppolys().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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)  ppoly_E,
real, dimension(:,:), intent(inout)  ppoly_S,
real, dimension(:,:), intent(inout)  ppoly_coef,
real, intent(in), optional  h_neglect 
)
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]ppoly_eEdge 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]

Definition at line 61 of file P3M_functions.F90.

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 

References regrid_edge_values::average_discontinuous_edge_values(), regrid_edge_values::bound_edge_values(), build_cubic_interpolant(), hneglect_dflt, is_cubic_monotonic(), and monotonize_cubic().

Referenced by p3m_interpolation().

Here is the call graph for this function:
Here is the caller graph for this function:

Variable Documentation

◆ hneglect_dflt

real, parameter p3m_functions::hneglect_dflt = 1.E-30
private

Default value of a negligible cell thickness.

Definition at line 13 of file P3M_functions.F90.

13 real, parameter :: hNeglect_dflt = 1.e-30 !< Default value of a negligible cell thickness

Referenced by p3m_boundary_extrapolation(), and p3m_limiter().

◆ hneglect_edge_dflt

real, parameter p3m_functions::hneglect_edge_dflt = 1.E-10
private

Default value of a negligible edge thickness.

Definition at line 14 of file P3M_functions.F90.

14 real, parameter :: hNeglect_edge_dflt = 1.e-10 !< Default value of a negligible edge thickness

Referenced by p3m_boundary_extrapolation().