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, edge_values, edge_slopes, ppoly_coef, h_neglect, answers_2018)
 Reconstruction by quartic polynomials within each cell. More...
 
subroutine pqm_limiter (N, h, u, edge_values, edge_slopes, h_neglect, answers_2018)
 Limit the piecewise quartic method reconstruction. More...
 
subroutine, public pqm_boundary_extrapolation (N, h, u, edge_values, ppoly_coef)
 Reconstruction by parabolas within boundary cells. More...
 
subroutine, public pqm_boundary_extrapolation_v1 (N, h, u, edge_values, edge_slopes, ppoly_coef, h_neglect)
 Reconstruction by parabolas within boundary cells. More...
 

Variables

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

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)  edge_values,
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) [H]
[in]ucell averages (size N) [A]
[in,out]edge_valuesEdge value of polynomial [A]
[in,out]ppoly_coefCoefficients of polynomial, mainly [A]

Definition at line 354 of file PQM_functions.F90.

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

◆ 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)  edge_values,
real, dimension(:,:), intent(inout)  edge_slopes,
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) [H]
[in]ucell averages (size N) [A]
[in,out]edge_valuesEdge value of polynomial [A]
[in,out]edge_slopesEdge slope of polynomial [A H-1]
[in,out]ppoly_coefCoefficients of polynomial, mainly [A]
[in]h_neglectA negligibly small width for the purpose of cell reconstructions [H]

Definition at line 501 of file PQM_functions.F90.

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

◆ 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)  edge_values,
real, dimension(:,:), intent(inout)  edge_slopes,
real, intent(in), optional  h_neglect,
logical, intent(in), optional  answers_2018 
)
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) [H]
[in]ucell average properties (size N) [A]
[in,out]edge_valuesPotentially modified edge values [A]
[in,out]edge_slopesPotentially modified edge slopes [A H-1]
[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 75 of file PQM_functions.F90.

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

◆ 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)  edge_values,
real, dimension(:,:), intent(inout)  edge_slopes,
real, dimension(:,:), intent(inout)  ppoly_coef,
real, intent(in), optional  h_neglect,
logical, intent(in), optional  answers_2018 
)

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) [H]
[in]ucell averages (size N) [A]
[in,out]edge_valuesEdge value of polynomial [A]
[in,out]edge_slopesEdge slope of polynomial [A H-1]
[in,out]ppoly_coefCoefficients of polynomial, mainly [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 20 of file PQM_functions.F90.

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