MOM6
Functions/Subroutines | Variables
pqm_functions Module Reference

Detailed Description

Piecewise quartic reconstruction functions.

Date of creation: 2008.06.06 L. White

This module contains routines that handle one-dimensionnal finite volume reconstruction using the piecewise quartic method (PQM).

Functions/Subroutines

subroutine, public pqm_reconstruction (N, h, u, ppoly_E, ppoly_S, ppoly_coef, h_neglect)
 Reconstruction by quartic polynomials within each cell. More...
 
subroutine pqm_limiter (N, h, u, ppoly_E, ppoly_S, h_neglect)
 Limit the piecewise quartic method reconstruction. More...
 
subroutine, public pqm_boundary_extrapolation (N, h, u, ppoly_E, ppoly_coef)
 Reconstruction by parabolas within boundary cells. More...
 
subroutine, public pqm_boundary_extrapolation_v1 (N, h, u, ppoly_E, ppoly_S, ppoly_coef, h_neglect)
 Reconstruction by parabolas within boundary cells. More...
 

Variables

real, parameter hneglect_dflt = 1.E-30
 Default negligible cell thickness. More...
 

Function/Subroutine Documentation

◆ pqm_boundary_extrapolation()

subroutine, public pqm_functions::pqm_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_coef 
)

Reconstruction by parabolas within boundary cells.

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

A parabola needs to be built in the cell and requires three degrees of freedom, which are the right edge value and slope and the cell average. 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 resulting parabola is not necessarily monotonic and the traditional PPM limiter is used to modify one of the edge values in order to yield a monotonic parabola.

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)
[in]ucell averages (size N)
[in,out]ppoly_eEdge value of polynomial, with the same units as u.
[in,out]ppoly_coefCoefficients of polynomial, mainly with the same units as u.

Definition at line 359 of file PQM_functions.F90.

359  integer, intent(in) :: N !< Number of cells
360  real, dimension(:), intent(in) :: h !< cell widths (size N)
361  real, dimension(:), intent(in) :: u !< cell averages (size N)
362  real, dimension(:,:), intent(inout) :: ppoly_E !< Edge value of polynomial,
363  !! with the same units as u.
364  real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial, mainly
365  !! with the same units as u.
366  ! Local variables
367  integer :: i0, i1
368  real :: u0, u1
369  real :: h0, h1
370  real :: a, b, c, d, e
371  real :: u0_l, u0_r
372  real :: u1_l, u1_r
373  real :: slope
374  real :: exp1, exp2
375 
376  ! ----- Left boundary -----
377  i0 = 1
378  i1 = 2
379  h0 = h(i0)
380  h1 = h(i1)
381  u0 = u(i0)
382  u1 = u(i1)
383 
384  ! Compute the left edge slope in neighboring cell and express it in
385  ! the global coordinate system
386  b = ppoly_coef(i1,2)
387  u1_r = b *(h0/h1) ! derivative evaluated at xi = 0.0,
388  ! expressed w.r.t. xi (local coord. system)
389 
390  ! Limit the right slope by the PLM limited slope
391  slope = 2.0 * ( u1 - u0 )
392  if ( abs(u1_r) > abs(slope) ) then
393  u1_r = slope
394  endif
395 
396  ! The right edge value in the boundary cell is taken to be the left
397  ! edge value in the neighboring cell
398  u0_r = ppoly_e(i1,1)
399 
400  ! Given the right edge value and slope, we determine the left
401  ! edge value and slope by computing the parabola as determined by
402  ! the right edge value and slope and the boundary cell average
403  u0_l = 3.0 * u0 + 0.5 * u1_r - 2.0 * u0_r
404 
405  ! Apply the traditional PPM limiter
406  exp1 = (u0_r - u0_l) * (u0 - 0.5*(u0_l+u0_r))
407  exp2 = (u0_r - u0_l) * (u0_r - u0_l) / 6.0
408 
409  if ( exp1 > exp2 ) then
410  u0_l = 3.0 * u0 - 2.0 * u0_r
411  endif
412 
413  if ( exp1 < -exp2 ) then
414  u0_r = 3.0 * u0 - 2.0 * u0_l
415  endif
416 
417  ppoly_e(i0,1) = u0_l
418  ppoly_e(i0,2) = u0_r
419 
420  a = u0_l
421  b = 6.0 * u0 - 4.0 * u0_l - 2.0 * u0_r
422  c = 3.0 * ( u0_r + u0_l - 2.0 * u0 )
423 
424  ! The quartic is reduced to a parabola in the boundary cell
425  ppoly_coef(i0,1) = a
426  ppoly_coef(i0,2) = b
427  ppoly_coef(i0,3) = c
428  ppoly_coef(i0,4) = 0.0
429  ppoly_coef(i0,5) = 0.0
430 
431  ! ----- Right boundary -----
432  i0 = n-1
433  i1 = n
434  h0 = h(i0)
435  h1 = h(i1)
436  u0 = u(i0)
437  u1 = u(i1)
438 
439  ! Compute the right edge slope in neighboring cell and express it in
440  ! the global coordinate system
441  b = ppoly_coef(i0,2)
442  c = ppoly_coef(i0,3)
443  d = ppoly_coef(i0,4)
444  e = ppoly_coef(i0,5)
445  u1_l = (b + 2*c + 3*d + 4*e) ! derivative evaluated at xi = 1.0
446  u1_l = u1_l * (h1/h0)
447 
448  ! Limit the left slope by the PLM limited slope
449  slope = 2.0 * ( u1 - u0 )
450  if ( abs(u1_l) > abs(slope) ) then
451  u1_l = slope
452  endif
453 
454  ! The left edge value in the boundary cell is taken to be the right
455  ! edge value in the neighboring cell
456  u0_l = ppoly_e(i0,2)
457 
458  ! Given the left edge value and slope, we determine the right
459  ! edge value and slope by computing the parabola as determined by
460  ! the left edge value and slope and the boundary cell average
461  u0_r = 3.0 * u1 - 0.5 * u1_l - 2.0 * u0_l
462 
463  ! Apply the traditional PPM limiter
464  exp1 = (u0_r - u0_l) * (u1 - 0.5*(u0_l+u0_r))
465  exp2 = (u0_r - u0_l) * (u0_r - u0_l) / 6.0
466 
467  if ( exp1 > exp2 ) then
468  u0_l = 3.0 * u1 - 2.0 * u0_r
469  endif
470 
471  if ( exp1 < -exp2 ) then
472  u0_r = 3.0 * u1 - 2.0 * u0_l
473  endif
474 
475  ppoly_e(i1,1) = u0_l
476  ppoly_e(i1,2) = u0_r
477 
478  a = u0_l
479  b = 6.0 * u1 - 4.0 * u0_l - 2.0 * u0_r
480  c = 3.0 * ( u0_r + u0_l - 2.0 * u1 )
481 
482  ! The quartic is reduced to a parabola in the boundary cell
483  ppoly_coef(i1,1) = a
484  ppoly_coef(i1,2) = b
485  ppoly_coef(i1,3) = c
486  ppoly_coef(i1,4) = 0.0
487  ppoly_coef(i1,5) = 0.0
488 

◆ pqm_boundary_extrapolation_v1()

subroutine, public pqm_functions::pqm_boundary_extrapolation_v1 ( 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 
)

Reconstruction by parabolas within boundary cells.

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

A parabola needs to be built in the cell and requires three degrees of freedom, which are the right edge value and slope and the cell average. 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 resulting parabola is not necessarily monotonic and the traditional PPM limiter is used to modify one of the edge values in order to yield a monotonic parabola.

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)
[in]ucell averages (size N)
[in,out]ppoly_eEdge value of polynomial, with the same units as u.
[in,out]ppoly_sEdge slope of polynomial, in the units of u over the units of h.
[in,out]ppoly_coefCoefficients of polynomial, mainly with the same units as u.
[in]h_neglectA negligibly small width for the purpose of cell reconstructions in the same units as h.

Definition at line 508 of file PQM_functions.F90.

508  integer, intent(in) :: N !< Number of cells
509  real, dimension(:), intent(in) :: h !< cell widths (size N)
510  real, dimension(:), intent(in) :: u !< cell averages (size N)
511  real, dimension(:,:), intent(inout) :: ppoly_E !< Edge value of polynomial,
512  !! with the same units as u.
513  real, dimension(:,:), intent(inout) :: ppoly_S !< Edge slope of polynomial,
514  !! in the units of u over the units of h.
515  real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial, mainly
516  !! with the same units as u.
517  real, optional, intent(in) :: h_neglect !< A negligibly small width for
518  !! the purpose of cell reconstructions
519  !! in the same units as h.
520  ! Local variables
521  integer :: i0, i1
522  integer :: inflexion_l
523  integer :: inflexion_r
524  real :: u0, u1, um
525  real :: h0, h1
526  real :: a, b, c, d, e
527  real :: ar, br, beta
528  real :: u0_l, u0_r
529  real :: u1_l, u1_r
530  real :: u_plm
531  real :: slope
532  real :: alpha1, alpha2, alpha3
533  real :: rho, sqrt_rho
534  real :: gradient1, gradient2
535  real :: x1, x2
536  real :: hNeglect
537 
538  hneglect = hneglect_dflt ; if (present(h_neglect)) hneglect = h_neglect
539 
540  ! ----- Left boundary (TOP) -----
541  i0 = 1
542  i1 = 2
543  h0 = h(i0)
544  h1 = h(i1)
545  u0 = u(i0)
546  u1 = u(i1)
547  um = u0
548 
549  ! Compute real slope and express it w.r.t. local coordinate system
550  ! within boundary cell
551  slope = 2.0 * ( u1 - u0 ) / ( ( h0 + h1 ) + hneglect )
552  slope = slope * h0
553 
554  ! The right edge value and slope of the boundary cell are taken to be the
555  ! left edge value and slope of the adjacent cell
556  a = ppoly_coef(i1,1)
557  b = ppoly_coef(i1,2)
558 
559  u0_r = a ! edge value
560  u1_r = b / (h1 + hneglect) ! edge slope (w.r.t. global coord.)
561 
562  ! Compute coefficient for rational function based on mean and right
563  ! edge value and slope
564  if (u1_r /= 0.) then ! HACK by AJA
565  beta = 2.0 * ( u0_r - um ) / ( (h0 + hneglect)*u1_r) - 1.0
566  else
567  beta = 0.
568  endif ! HACK by AJA
569  br = u0_r + beta*u0_r - um
570  ar = um + beta*um - br
571 
572  ! Left edge value estimate based on rational function
573  u0_l = ar
574 
575  ! Edge value estimate based on PLM
576  u_plm = um - 0.5 * slope
577 
578  ! Check whether the left edge value is bounded by the mean and
579  ! the PLM edge value. If so, keep it and compute left edge slope
580  ! based on the rational function. If not, keep the PLM edge value and
581  ! compute corresponding slope.
582  if ( abs(um-u0_l) < abs(um-u_plm) ) then
583  u1_l = 2.0 * ( br - ar*beta)
584  u1_l = u1_l / (h0 + hneglect)
585  else
586  u0_l = u_plm
587  u1_l = slope / (h0 + hneglect)
588  endif
589 
590  ! Monotonize quartic
591  inflexion_l = 0
592 
593  a = u0_l
594  b = h0 * u1_l
595  c = 30.0 * um - 12.0*u0_r - 18.0*u0_l + 1.5*h0*(u1_r - 3.0*u1_l)
596  d = -60.0 * um + h0 *(6.0*u1_l - 4.0*u1_r) + 28.0*u0_r + 32.0*u0_l
597  e = 30.0 * um + 2.5*h0*(u1_r - u1_l) - 15.0*(u0_l + u0_r)
598 
599  alpha1 = 6*e
600  alpha2 = 3*d
601  alpha3 = c
602 
603  rho = alpha2 * alpha2 - 4.0 * alpha1 * alpha3
604 
605  ! Check whether inflexion points exist. If so, transform the quartic
606  ! so that both inflexion points coalesce on the left edge.
607  if (( alpha1 /= 0.0 ) .and. ( rho >= 0.0 )) then
608 
609  sqrt_rho = sqrt( rho )
610 
611  x1 = 0.5 * ( - alpha2 - sqrt_rho ) / alpha1
612  if ( (x1 > 0.0) .and. (x1 < 1.0) ) then
613  gradient1 = 4.0 * e * (x1**3) + 3.0 * d * (x1**2) + 2.0 * c * x1 + b
614  if ( gradient1 * slope < 0.0 ) then
615  inflexion_l = 1
616  endif
617  endif
618 
619  x2 = 0.5 * ( - alpha2 + sqrt_rho ) / alpha1
620  if ( (x2 > 0.0) .and. (x2 < 1.0) ) then
621  gradient2 = 4.0 * e * (x2**3) + 3.0 * d * (x2**2) + 2.0 * c * x2 + b
622  if ( gradient2 * slope < 0.0 ) then
623  inflexion_l = 1
624  endif
625  endif
626 
627  endif
628 
629  if (( alpha1 == 0.0 ) .and. ( alpha2 /= 0.0 )) then
630 
631  x1 = - alpha3 / alpha2
632  if ( (x1 >= 0.0) .and. (x1 <= 1.0) ) then
633  gradient1 = 3.0 * d * (x1**2) + 2.0 * c * x1 + b
634  if ( gradient1 * slope < 0.0 ) then
635  inflexion_l = 1
636  endif
637  endif
638 
639  endif
640 
641  if ( inflexion_l == 1 ) then
642 
643  ! We modify the edge slopes so that both inflexion points
644  ! collapse onto the left edge
645  u1_l = ( 10.0 * um - 2.0 * u0_r - 8.0 * u0_l ) / (3.0*h0 + hneglect)
646  u1_r = ( -10.0 * um + 6.0 * u0_r + 4.0 * u0_l ) / (h0 + hneglect)
647 
648  ! One of the modified slopes might be inconsistent. When that happens,
649  ! the inconsistent slope is set equal to zero and the opposite edge value
650  ! and edge slope are modified in compliance with the fact that both
651  ! inflexion points must still be located on the left edge
652  if ( u1_l * slope < 0.0 ) then
653 
654  u1_l = 0.0
655  u0_r = 5.0 * um - 4.0 * u0_l
656  u1_r = 20.0 * (um - u0_l) / ( h0 + hneglect )
657 
658  elseif ( u1_r * slope < 0.0 ) then
659 
660  u1_r = 0.0
661  u0_l = (5.0*um - 3.0*u0_r) / 2.0
662  u1_l = 10.0 * (-um + u0_r) / (3.0 * h0 + hneglect )
663 
664  endif
665 
666  endif
667 
668  ! Store edge values, edge slopes and coefficients
669  ppoly_e(i0,1) = u0_l
670  ppoly_e(i0,2) = u0_r
671  ppoly_s(i0,1) = u1_l
672  ppoly_s(i0,2) = u1_r
673 
674  a = u0_l
675  b = h0 * u1_l
676  c = 30.0 * um - 12.0*u0_r - 18.0*u0_l + 1.5*h0*(u1_r - 3.0*u1_l)
677  d = -60.0 * um + h0 *(6.0*u1_l - 4.0*u1_r) + 28.0*u0_r + 32.0*u0_l
678  e = 30.0 * um + 2.5*h0*(u1_r - u1_l) - 15.0*(u0_l + u0_r)
679 
680  ! Store coefficients
681  ppoly_coef(i0,1) = a
682  ppoly_coef(i0,2) = b
683  ppoly_coef(i0,3) = c
684  ppoly_coef(i0,4) = d
685  ppoly_coef(i0,5) = e
686 
687  ! ----- Right boundary (BOTTOM) -----
688  i0 = n-1
689  i1 = n
690  h0 = h(i0)
691  h1 = h(i1)
692  u0 = u(i0)
693  u1 = u(i1)
694  um = u1
695 
696  ! Compute real slope and express it w.r.t. local coordinate system
697  ! within boundary cell
698  slope = 2.0 * ( u1 - u0 ) / ( h0 + h1 )
699  slope = slope * h1
700 
701  ! The left edge value and slope of the boundary cell are taken to be the
702  ! right edge value and slope of the adjacent cell
703  a = ppoly_coef(i0,1)
704  b = ppoly_coef(i0,2)
705  c = ppoly_coef(i0,3)
706  d = ppoly_coef(i0,4)
707  e = ppoly_coef(i0,5)
708  u0_l = a + b + c + d + e ! edge value
709  u1_l = (b + 2*c + 3*d + 4*e) / h0 ! edge slope (w.r.t. global coord.)
710 
711  ! Compute coefficient for rational function based on mean and left
712  ! edge value and slope
713  if (um-u0_l /= 0.) then ! HACK by AJA
714  beta = 0.5*h1*u1_l / (um-u0_l) - 1.0
715  else
716  beta = 0.
717  endif ! HACK by AJA
718  br = beta*um + um - u0_l
719  ar = u0_l
720 
721  ! Right edge value estimate based on rational function
722  if (1+beta /= 0.) then ! HACK by AJA
723  u0_r = (ar + 2*br + beta*br ) / ((1+beta)*(1+beta))
724  else
725  u0_r = um + 0.5 * slope ! PLM
726  endif ! HACK by AJA
727 
728  ! Right edge value estimate based on PLM
729  u_plm = um + 0.5 * slope
730 
731  ! Check whether the right edge value is bounded by the mean and
732  ! the PLM edge value. If so, keep it and compute right edge slope
733  ! based on the rational function. If not, keep the PLM edge value and
734  ! compute corresponding slope.
735  if ( abs(um-u0_r) < abs(um-u_plm) ) then
736  u1_r = 2.0 * ( br - ar*beta ) / ( (1+beta)*(1+beta)*(1+beta) )
737  u1_r = u1_r / h1
738  else
739  u0_r = u_plm
740  u1_r = slope / h1
741  endif
742 
743  ! Monotonize quartic
744  inflexion_r = 0
745 
746  a = u0_l
747  b = h1 * u1_l
748  c = 30.0 * um - 12.0*u0_r - 18.0*u0_l + 1.5*h1*(u1_r - 3.0*u1_l)
749  d = -60.0 * um + h1*(6.0*u1_l - 4.0*u1_r) + 28.0*u0_r + 32.0*u0_l
750  e = 30.0 * um + 2.5*h1*(u1_r - u1_l) - 15.0*(u0_l + u0_r)
751 
752  alpha1 = 6*e
753  alpha2 = 3*d
754  alpha3 = c
755 
756  rho = alpha2 * alpha2 - 4.0 * alpha1 * alpha3
757 
758  ! Check whether inflexion points exist. If so, transform the quartic
759  ! so that both inflexion points coalesce on the right edge.
760  if (( alpha1 /= 0.0 ) .and. ( rho >= 0.0 )) then
761 
762  sqrt_rho = sqrt( rho )
763 
764  x1 = 0.5 * ( - alpha2 - sqrt_rho ) / alpha1
765  if ( (x1 > 0.0) .and. (x1 < 1.0) ) then
766  gradient1 = 4.0 * e * (x1**3) + 3.0 * d * (x1**2) + 2.0 * c * x1 + b
767  if ( gradient1 * slope < 0.0 ) then
768  inflexion_r = 1
769  endif
770  endif
771 
772  x2 = 0.5 * ( - alpha2 + sqrt_rho ) / alpha1
773  if ( (x2 > 0.0) .and. (x2 < 1.0) ) then
774  gradient2 = 4.0 * e * (x2**3) + 3.0 * d * (x2**2) + 2.0 * c * x2 + b
775  if ( gradient2 * slope < 0.0 ) then
776  inflexion_r = 1
777  endif
778  endif
779 
780  endif
781 
782  if (( alpha1 == 0.0 ) .and. ( alpha2 /= 0.0 )) then
783 
784  x1 = - alpha3 / alpha2
785  if ( (x1 >= 0.0) .and. (x1 <= 1.0) ) then
786  gradient1 = 3.0 * d * (x1**2) + 2.0 * c * x1 + b
787  if ( gradient1 * slope < 0.0 ) then
788  inflexion_r = 1
789  endif
790  endif
791 
792  endif
793 
794  if ( inflexion_r == 1 ) then
795 
796  ! We modify the edge slopes so that both inflexion points
797  ! collapse onto the right edge
798  u1_r = ( -10.0 * um + 8.0 * u0_r + 2.0 * u0_l ) / (3.0 * h1)
799  u1_l = ( 10.0 * um - 4.0 * u0_r - 6.0 * u0_l ) / h1
800 
801  ! One of the modified slopes might be inconsistent. When that happens,
802  ! the inconsistent slope is set equal to zero and the opposite edge value
803  ! and edge slope are modified in compliance with the fact that both
804  ! inflexion points must still be located on the right edge
805  if ( u1_l * slope < 0.0 ) then
806 
807  u1_l = 0.0
808  u0_r = ( 5.0 * um - 3.0 * u0_l ) / 2.0
809  u1_r = 10.0 * (um - u0_l) / (3.0 * h1)
810 
811  elseif ( u1_r * slope < 0.0 ) then
812 
813  u1_r = 0.0
814  u0_l = 5.0 * um - 4.0 * u0_r
815  u1_l = 20.0 * ( -um + u0_r ) / h1
816 
817  endif
818 
819  endif
820 
821  ! Store edge values, edge slopes and coefficients
822  ppoly_e(i1,1) = u0_l
823  ppoly_e(i1,2) = u0_r
824  ppoly_s(i1,1) = u1_l
825  ppoly_s(i1,2) = u1_r
826 
827  a = u0_l
828  b = h1 * u1_l
829  c = 30.0 * um - 12.0*u0_r - 18.0*u0_l + 1.5*h1*(u1_r - 3.0*u1_l)
830  d = -60.0 * um + h1 *(6.0*u1_l - 4.0*u1_r) + 28.0*u0_r + 32.0*u0_l
831  e = 30.0 * um + 2.5*h1*(u1_r - u1_l) - 15.0*(u0_l + u0_r)
832 
833  ppoly_coef(i1,1) = a
834  ppoly_coef(i1,2) = b
835  ppoly_coef(i1,3) = c
836  ppoly_coef(i1,4) = d
837  ppoly_coef(i1,5) = e
838 

References hneglect_dflt.

Referenced by mom_remapping::build_reconstructions_1d(), and regrid_interp::regridding_set_ppolys().

Here is the caller graph for this function:

◆ pqm_limiter()

subroutine pqm_functions::pqm_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, intent(in), optional  h_neglect 
)
private

Limit the piecewise quartic method reconstruction.

Standard PQM limiter (White & Adcroft, JCP 2008).

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

Parameters
[in]nNumber of cells
[in]hcell widths (size N)
[in]ucell average properties (size N)
[in,out]ppoly_ePotentially modified edge values, with the same units as u.
[in,out]ppoly_sPotentially modified edge slopes, with the same units as u.
[in]h_neglectA negligibly small width for the purpose of cell reconstructions in the same units as h

Definition at line 79 of file PQM_functions.F90.

79  integer, intent(in) :: N !< Number of cells
80  real, dimension(:), intent(in) :: h !< cell widths (size N)
81  real, dimension(:), intent(in) :: u !< cell average properties (size N)
82  real, dimension(:,:), intent(inout) :: ppoly_E !< Potentially modified edge values,
83  !! with the same units as u.
84  real, dimension(:,:), intent(inout) :: ppoly_S !< Potentially modified edge slopes,
85  !! with the same units as u.
86  real, optional, intent(in) :: h_neglect !< A negligibly small width for
87  !! the purpose of cell reconstructions
88  !! in the same units as h
89  ! Local variables
90  integer :: k ! loop index
91  integer :: inflexion_l
92  integer :: inflexion_r
93  real :: u0_l, u0_r ! edge values
94  real :: u1_l, u1_r ! edge slopes
95  real :: u_l, u_c, u_r ! left, center and right cell averages
96  real :: h_l, h_c, h_r ! left, center and right cell widths
97  real :: sigma_l, sigma_c, sigma_r ! left, center and right van Leer slopes
98  real :: slope ! retained PLM slope
99  real :: a, b, c, d, e
100  real :: alpha1, alpha2, alpha3
101  real :: rho, sqrt_rho
102  real :: gradient1, gradient2
103  real :: x1, x2
104  real :: hNeglect
105 
106  hneglect = hneglect_dflt ; if (present(h_neglect)) hneglect = h_neglect
107 
108  ! Bound edge values
109  call bound_edge_values( n, h, u, ppoly_e, hneglect )
110 
111  ! Make discontinuous edge values monotonic (thru averaging)
112  call check_discontinuous_edge_values( n, u, ppoly_e )
113 
114  ! Loop on interior cells to apply the PQM limiter
115  do k = 2,n-1
116 
117  !if ( h(k) < 1.0 ) cycle
118 
119  inflexion_l = 0
120  inflexion_r = 0
121 
122  ! Get edge values, edge slopes and cell width
123  u0_l = ppoly_e(k,1)
124  u0_r = ppoly_e(k,2)
125  u1_l = ppoly_s(k,1)
126  u1_r = ppoly_s(k,2)
127 
128  ! Get cell widths and cell averages (boundary cells are assumed to
129  ! be local extrema for the sake of slopes)
130  h_l = h(k-1)
131  h_c = h(k)
132  h_r = h(k+1)
133  u_l = u(k-1)
134  u_c = u(k)
135  u_r = u(k+1)
136 
137  ! Compute limited slope
138  sigma_l = 2.0 * ( u_c - u_l ) / ( h_c + hneglect )
139  sigma_c = 2.0 * ( u_r - u_l ) / ( h_l + 2.0*h_c + h_r + hneglect )
140  sigma_r = 2.0 * ( u_r - u_c ) / ( h_c + hneglect )
141 
142  if ( (sigma_l * sigma_r) > 0.0 ) then
143  slope = sign( min(abs(sigma_l),abs(sigma_c),abs(sigma_r)), sigma_c )
144  else
145  slope = 0.0
146  endif
147 
148  ! If one of the slopes has the wrong sign compared with the
149  ! limited PLM slope, it is set equal to the limited PLM slope
150  if ( u1_l*slope <= 0.0 ) u1_l = slope
151  if ( u1_r*slope <= 0.0 ) u1_r = slope
152 
153  ! Local extremum --> flatten
154  if ( (u0_r - u_c) * (u_c - u0_l) <= 0.0) then
155  u0_l = u_c
156  u0_r = u_c
157  u1_l = 0.0
158  u1_r = 0.0
159  inflexion_l = -1
160  inflexion_r = -1
161  endif
162 
163  ! Edge values are bounded and averaged when discontinuous and not
164  ! monotonic, edge slopes are consistent and the cell is not an extremum.
165  ! We now need to check and encorce the monotonicity of the quartic within
166  ! the cell
167  if ( (inflexion_l == 0) .AND. (inflexion_r == 0) ) then
168 
169  a = u0_l
170  b = h_c * u1_l
171  c = 30.0 * u(k) - 12.0*u0_r - 18.0*u0_l + 1.5*h_c*(u1_r - 3.0*u1_l)
172  d = -60.0 * u(k) + h_c *(6.0*u1_l - 4.0*u1_r) + 28.0*u0_r + 32.0*u0_l
173  e = 30.0 * u(k) + 2.5*h_c*(u1_r - u1_l) - 15.0*(u0_l + u0_r)
174 
175  ! Determine the coefficients of the second derivative
176  ! alpha1 xi^2 + alpha2 xi + alpha3
177  alpha1 = 6*e
178  alpha2 = 3*d
179  alpha3 = c
180 
181  rho = alpha2 * alpha2 - 4.0 * alpha1 * alpha3
182 
183  ! Check whether inflexion points exist
184  if (( alpha1 /= 0.0 ) .and. ( rho >= 0.0 )) then
185 
186  sqrt_rho = sqrt( rho )
187 
188  x1 = 0.5 * ( - alpha2 - sqrt_rho ) / alpha1
189  x2 = 0.5 * ( - alpha2 + sqrt_rho ) / alpha1
190 
191  ! Check whether both inflexion points lie in [0,1]
192  if ( (x1 >= 0.0) .AND. (x1 <= 1.0) .AND. &
193  (x2 >= 0.0) .AND. (x2 <= 1.0) ) then
194 
195  gradient1 = 4.0 * e * (x1**3) + 3.0 * d * (x1**2) + 2.0 * c * x1 + b
196  gradient2 = 4.0 * e * (x2**3) + 3.0 * d * (x2**2) + 2.0 * c * x2 + b
197 
198  ! Check whether one of the gradients is inconsistent
199  if ( (gradient1 * slope < 0.0) .OR. &
200  (gradient2 * slope < 0.0) ) then
201  ! Decide where to collapse inflexion points
202  ! (depends on one-sided slopes)
203  if ( abs(sigma_l) < abs(sigma_r) ) then
204  inflexion_l = 1
205  else
206  inflexion_r = 1
207  endif
208  endif
209 
210  ! If both x1 and x2 do not lie in [0,1], check whether
211  ! only x1 lies in [0,1]
212  elseif ( (x1 >= 0.0) .AND. (x1 <= 1.0) ) then
213 
214  gradient1 = 4.0 * e * (x1**3) + 3.0 * d * (x1**2) + 2.0 * c * x1 + b
215 
216  ! Check whether the gradient is inconsistent
217  if ( gradient1 * slope < 0.0 ) then
218  ! Decide where to collapse inflexion points
219  ! (depends on one-sided slopes)
220  if ( abs(sigma_l) < abs(sigma_r) ) then
221  inflexion_l = 1
222  else
223  inflexion_r = 1
224  endif
225  endif
226 
227  ! If x1 does not lie in [0,1], check whether x2 lies in [0,1]
228  elseif ( (x2 >= 0.0) .AND. (x2 <= 1.0) ) then
229 
230  gradient2 = 4.0 * e * (x2**3) + 3.0 * d * (x2**2) + 2.0 * c * x2 + b
231 
232  ! Check whether the gradient is inconsistent
233  if ( gradient2 * slope < 0.0 ) then
234  ! Decide where to collapse inflexion points
235  ! (depends on one-sided slopes)
236  if ( abs(sigma_l) < abs(sigma_r) ) then
237  inflexion_l = 1
238  else
239  inflexion_r = 1
240  endif
241  endif
242 
243  endif ! end checking where the inflexion points lie
244 
245  endif ! end checking if alpha1 != 0 AND rho >= 0
246 
247  ! If alpha1 is zero, the second derivative of the quartic reduces
248  ! to a straight line
249  if (( alpha1 == 0.0 ) .and. ( alpha2 /= 0.0 )) then
250 
251  x1 = - alpha3 / alpha2
252  if ( (x1 >= 0.0) .AND. (x1 <= 1.0) ) then
253 
254  gradient1 = 4.0 * e * (x1**3) + 3.0 * d * (x1**2) + 2.0 * c * x1 + b
255 
256  ! Check whether the gradient is inconsistent
257  if ( gradient1 * slope < 0.0 ) then
258  ! Decide where to collapse inflexion points
259  ! (depends on one-sided slopes)
260  if ( abs(sigma_l) < abs(sigma_r) ) then
261  inflexion_l = 1
262  else
263  inflexion_r = 1
264  endif
265  endif ! check slope consistency
266 
267  endif
268 
269  endif ! end check whether we can find the root of the straight line
270 
271  endif ! end checking whether to shift inflexion points
272 
273  ! At this point, we know onto which edge to shift inflexion points
274  if ( inflexion_l == 1 ) then
275 
276  ! We modify the edge slopes so that both inflexion points
277  ! collapse onto the left edge
278  u1_l = ( 10.0 * u_c - 2.0 * u0_r - 8.0 * u0_l ) / (3.0*h_c + hneglect )
279  u1_r = ( -10.0 * u_c + 6.0 * u0_r + 4.0 * u0_l ) / ( h_c + hneglect )
280 
281  ! One of the modified slopes might be inconsistent. When that happens,
282  ! the inconsistent slope is set equal to zero and the opposite edge value
283  ! and edge slope are modified in compliance with the fact that both
284  ! inflexion points must still be located on the left edge
285  if ( u1_l * slope < 0.0 ) then
286 
287  u1_l = 0.0
288  u0_r = 5.0 * u_c - 4.0 * u0_l
289  u1_r = 20.0 * (u_c - u0_l) / ( h_c + hneglect )
290 
291  elseif ( u1_r * slope < 0.0 ) then
292 
293  u1_r = 0.0
294  u0_l = (5.0*u_c - 3.0*u0_r) / 2.0
295  u1_l = 10.0 * (-u_c + u0_r) / (3.0 * h_c + hneglect)
296 
297  endif
298 
299  elseif ( inflexion_r == 1 ) then
300 
301  ! We modify the edge slopes so that both inflexion points
302  ! collapse onto the right edge
303  u1_r = ( -10.0 * u_c + 8.0 * u0_r + 2.0 * u0_l ) / (3.0 * h_c + hneglect)
304  u1_l = ( 10.0 * u_c - 4.0 * u0_r - 6.0 * u0_l ) / (h_c + hneglect)
305 
306  ! One of the modified slopes might be inconsistent. When that happens,
307  ! the inconsistent slope is set equal to zero and the opposite edge value
308  ! and edge slope are modified in compliance with the fact that both
309  ! inflexion points must still be located on the right edge
310  if ( u1_l * slope < 0.0 ) then
311 
312  u1_l = 0.0
313  u0_r = ( 5.0 * u_c - 3.0 * u0_l ) / 2.0
314  u1_r = 10.0 * (u_c - u0_l) / (3.0 * h_c + hneglect)
315 
316  elseif ( u1_r * slope < 0.0 ) then
317 
318  u1_r = 0.0
319  u0_l = 5.0 * u_c - 4.0 * u0_r
320  u1_l = 20.0 * ( -u_c + u0_r ) / (h_c + hneglect)
321 
322  endif
323 
324  endif ! clause to check where to collapse inflexion points
325 
326  ! Save edge values and edge slopes for reconstruction
327  ppoly_e(k,1) = u0_l
328  ppoly_e(k,2) = u0_r
329  ppoly_s(k,1) = u1_l
330  ppoly_s(k,2) = u1_r
331 
332  enddo ! end loop on interior cells
333 
334  ! Constant reconstruction within boundary cells
335  ppoly_e(1,:) = u(1)
336  ppoly_s(1,:) = 0.0
337 
338  ppoly_e(n,:) = u(n)
339  ppoly_s(n,:) = 0.0
340 

References regrid_edge_values::bound_edge_values(), regrid_edge_values::check_discontinuous_edge_values(), and hneglect_dflt.

Referenced by pqm_reconstruction().

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

◆ pqm_reconstruction()

subroutine, public pqm_functions::pqm_reconstruction ( 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 
)

Reconstruction by quartic polynomials within each cell.

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

Parameters
[in]nNumber of cells
[in]hcell widths (size N)
[in]ucell averages (size N)
[in,out]ppoly_eEdge value of polynomial, with the same units as u.
[in,out]ppoly_sEdge slope of polynomial, in the units of u over the units of h.
[in,out]ppoly_coefCoefficients of polynomial, mainly with the same units as u.
[in]h_neglectA negligibly small width for the purpose of cell reconstructions in the same units as h

Definition at line 21 of file PQM_functions.F90.

21  integer, intent(in) :: N !< Number of cells
22  real, dimension(:), intent(in) :: h !< cell widths (size N)
23  real, dimension(:), intent(in) :: u !< cell averages (size N)
24  real, dimension(:,:), intent(inout) :: ppoly_E !< Edge value of polynomial,
25  !! with the same units as u.
26  real, dimension(:,:), intent(inout) :: ppoly_S !< Edge slope of polynomial,
27  !! in the units of u over the units of h.
28  real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial, mainly
29  !! with the same units as u.
30  real, optional, intent(in) :: h_neglect !< A negligibly small width for
31  !! the purpose of cell reconstructions
32  !! in the same units as h
33 
34  ! Local variables
35  integer :: k ! loop index
36  real :: h_c ! cell width
37  real :: u0_l, u0_r ! edge values (left and right)
38  real :: u1_l, u1_r ! edge slopes (left and right)
39  real :: a, b, c, d, e ! parabola coefficients
40 
41  ! PQM limiter
42  call pqm_limiter( n, h, u, ppoly_e, ppoly_s, h_neglect )
43 
44  ! Loop on cells to construct the cubic within each cell
45  do k = 1,n
46 
47  u0_l = ppoly_e(k,1)
48  u0_r = ppoly_e(k,2)
49 
50  u1_l = ppoly_s(k,1)
51  u1_r = ppoly_s(k,2)
52 
53  h_c = h(k)
54 
55  a = u0_l
56  b = h_c * u1_l
57  c = 30.0 * u(k) - 12.0*u0_r - 18.0*u0_l + 1.5*h_c*(u1_r - 3.0*u1_l)
58  d = -60.0 * u(k) + h_c *(6.0*u1_l - 4.0*u1_r) + 28.0*u0_r + 32.0*u0_l
59  e = 30.0 * u(k) + 2.5*h_c*(u1_r - u1_l) - 15.0*(u0_l + u0_r)
60 
61  ! Store coefficients
62  ppoly_coef(k,1) = a
63  ppoly_coef(k,2) = b
64  ppoly_coef(k,3) = c
65  ppoly_coef(k,4) = d
66  ppoly_coef(k,5) = e
67 
68  enddo ! end loop on cells
69 

References pqm_limiter().

Referenced by mom_remapping::build_reconstructions_1d(), and regrid_interp::regridding_set_ppolys().

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

Variable Documentation

◆ hneglect_dflt

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

Default negligible cell thickness.

Definition at line 12 of file PQM_functions.F90.

12 real, parameter :: hNeglect_dflt = 1.e-30 !< Default negligible cell thickness

Referenced by pqm_boundary_extrapolation_v1(), and pqm_limiter().