\hypertarget{namespacepqm__functions}{}\section{pqm\+\_\+functions Module Reference}
\label{namespacepqm__functions}\index{pqm\+\_\+functions@{pqm\+\_\+functions}}


\subsection{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 (P\+QM). \subsection*{Functions/\+Subroutines}
\begin{DoxyCompactItemize}
\item 
subroutine, public \mbox{\hyperlink{namespacepqm__functions_afa6f7b5430011f03c428b329c5f42fae}{pqm\+\_\+reconstruction}} (N, h, u, edge\+\_\+values, edge\+\_\+slopes, ppoly\+\_\+coef, h\+\_\+neglect, answers\+\_\+2018)
\begin{DoxyCompactList}\small\item\em Reconstruction by quartic polynomials within each cell. \end{DoxyCompactList}\item 
subroutine \mbox{\hyperlink{namespacepqm__functions_aab9f3108956943c39ed2d4b675d78021}{pqm\+\_\+limiter}} (N, h, u, edge\+\_\+values, edge\+\_\+slopes, h\+\_\+neglect, answers\+\_\+2018)
\begin{DoxyCompactList}\small\item\em Limit the piecewise quartic method reconstruction. \end{DoxyCompactList}\item 
subroutine, public \mbox{\hyperlink{namespacepqm__functions_ad970064d6c3540fd00dff2453d0778d7}{pqm\+\_\+boundary\+\_\+extrapolation}} (N, h, u, edge\+\_\+values, ppoly\+\_\+coef)
\begin{DoxyCompactList}\small\item\em Reconstruction by parabolas within boundary cells. \end{DoxyCompactList}\item 
subroutine, public \mbox{\hyperlink{namespacepqm__functions_a1e0e86d8470dd334b9cd676f511f6720}{pqm\+\_\+boundary\+\_\+extrapolation\+\_\+v1}} (N, h, u, edge\+\_\+values, edge\+\_\+slopes, ppoly\+\_\+coef, h\+\_\+neglect)
\begin{DoxyCompactList}\small\item\em Reconstruction by parabolas within boundary cells. \end{DoxyCompactList}\end{DoxyCompactItemize}
\subsection*{Variables}
\begin{DoxyCompactItemize}
\item 
\mbox{\Hypertarget{namespacepqm__functions_a3bd59442a660e0e66c29603c992932e6}\label{namespacepqm__functions_a3bd59442a660e0e66c29603c992932e6}} 
real, parameter \mbox{\hyperlink{namespacepqm__functions_a3bd59442a660e0e66c29603c992932e6}{hneglect\+\_\+dflt}} = 1.E-\/30
\begin{DoxyCompactList}\small\item\em Default negligible cell thickness. \end{DoxyCompactList}\end{DoxyCompactItemize}


\subsection{Function/\+Subroutine Documentation}
\mbox{\Hypertarget{namespacepqm__functions_ad970064d6c3540fd00dff2453d0778d7}\label{namespacepqm__functions_ad970064d6c3540fd00dff2453d0778d7}} 
\index{pqm\+\_\+functions@{pqm\+\_\+functions}!pqm\+\_\+boundary\+\_\+extrapolation@{pqm\+\_\+boundary\+\_\+extrapolation}}
\index{pqm\+\_\+boundary\+\_\+extrapolation@{pqm\+\_\+boundary\+\_\+extrapolation}!pqm\+\_\+functions@{pqm\+\_\+functions}}
\subsubsection{\texorpdfstring{pqm\+\_\+boundary\+\_\+extrapolation()}{pqm\_boundary\_extrapolation()}}
{\footnotesize\ttfamily subroutine, public pqm\+\_\+functions\+::pqm\+\_\+boundary\+\_\+extrapolation (\begin{DoxyParamCaption}\item[{integer, intent(in)}]{N,  }\item[{real, dimension(\+:), intent(in)}]{h,  }\item[{real, dimension(\+:), intent(in)}]{u,  }\item[{real, dimension(\+:,\+:), intent(inout)}]{edge\+\_\+values,  }\item[{real, dimension(\+:,\+:), intent(inout)}]{ppoly\+\_\+coef }\end{DoxyParamCaption})}



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 P\+PM 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 \textquotesingle{}u\textquotesingle{} is equal to the number of cells defining \textquotesingle{}grid\textquotesingle{} and \textquotesingle{}ppoly\textquotesingle{}. No consistency check is performed here.


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em n} & Number of cells\\
\hline
\mbox{\tt in}  & {\em h} & cell widths (size N) \mbox{[}H\mbox{]}\\
\hline
\mbox{\tt in}  & {\em u} & cell averages (size N) \mbox{[}A\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em edge\+\_\+values} & Edge value of polynomial \mbox{[}A\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em ppoly\+\_\+coef} & Coefficients of polynomial, mainly \mbox{[}A\mbox{]} \\
\hline
\end{DoxyParams}


Definition at line 355 of file P\+Q\+M\+\_\+functions.\+F90.


\begin{DoxyCode}
355   \textcolor{keywordtype}{integer},              \textcolor{keywordtype}{intent(in)}    :: N\textcolor{comment}{ !< Number of cells}
356   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:)},   \textcolor{keywordtype}{intent(in)}    :: h\textcolor{comment}{ !< cell widths (size N) [H]}
357   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:)},   \textcolor{keywordtype}{intent(in)}    :: u\textcolor{comment}{ !< cell averages (size N) [A]}
358   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:,:)}, \textcolor{keywordtype}{intent(inout)} :: edge\_values\textcolor{comment}{    !< Edge value of polynomial [A]}
359   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:,:)}, \textcolor{keywordtype}{intent(inout)} :: ppoly\_coef\textcolor{comment}{ !< Coefficients of polynomial, mainly [A]}
360   \textcolor{comment}{! Local variables}
361   \textcolor{keywordtype}{integer}       :: i0, i1
362   \textcolor{keywordtype}{real}          :: u0, u1
363   \textcolor{keywordtype}{real}          :: h0, h1
364   \textcolor{keywordtype}{real}          :: a, b, c, d, e
365   \textcolor{keywordtype}{real}          :: u0\_l, u0\_r
366   \textcolor{keywordtype}{real}          :: u1\_l, u1\_r
367   \textcolor{keywordtype}{real}          :: slope
368   \textcolor{keywordtype}{real}          :: exp1, exp2
369 
370   \textcolor{comment}{! ----- 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   \textcolor{comment}{! Compute the left edge slope in neighboring cell and express it in}
379   \textcolor{comment}{! the global coordinate system}
380   b = ppoly\_coef(i1,2)
381   u1\_r = b *(h0/h1)     \textcolor{comment}{! derivative evaluated at xi = 0.0,}
382                         \textcolor{comment}{! expressed w.r.t. xi (local coord. system)}
383 
384   \textcolor{comment}{! Limit the right slope by the PLM limited slope}
385   slope = 2.0 * ( u1 - u0 )
386   \textcolor{keywordflow}{if} ( abs(u1\_r) > abs(slope) ) \textcolor{keywordflow}{then}
387     u1\_r = slope
388 \textcolor{keywordflow}{  endif}
389 
390   \textcolor{comment}{! The right edge value in the boundary cell is taken to be the left}
391   \textcolor{comment}{! edge value in the neighboring cell}
392   u0\_r = edge\_values(i1,1)
393 
394   \textcolor{comment}{! Given the right edge value and slope, we determine the left}
395   \textcolor{comment}{! edge value and slope by computing the parabola as determined by}
396   \textcolor{comment}{! 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   \textcolor{comment}{! 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   \textcolor{keywordflow}{if} ( exp1 > exp2 ) \textcolor{keywordflow}{then}
404     u0\_l = 3.0 * u0 - 2.0 * u0\_r
405 \textcolor{keywordflow}{  endif}
406 
407   \textcolor{keywordflow}{if} ( exp1 < -exp2 ) \textcolor{keywordflow}{then}
408     u0\_r = 3.0 * u0 - 2.0 * u0\_l
409 \textcolor{keywordflow}{  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   \textcolor{comment}{! 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   \textcolor{comment}{! ----- 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   \textcolor{comment}{! Compute the right edge slope in neighboring cell and express it in}
434   \textcolor{comment}{! 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)      \textcolor{comment}{! derivative evaluated at xi = 1.0}
440   u1\_l = u1\_l * (h1/h0)
441 
442   \textcolor{comment}{! Limit the left slope by the PLM limited slope}
443   slope = 2.0 * ( u1 - u0 )
444   \textcolor{keywordflow}{if} ( abs(u1\_l) > abs(slope) ) \textcolor{keywordflow}{then}
445     u1\_l = slope
446 \textcolor{keywordflow}{  endif}
447 
448   \textcolor{comment}{! The left edge value in the boundary cell is taken to be the right}
449   \textcolor{comment}{! edge value in the neighboring cell}
450   u0\_l = edge\_values(i0,2)
451 
452   \textcolor{comment}{! Given the left edge value and slope, we determine the right}
453   \textcolor{comment}{! edge value and slope by computing the parabola as determined by}
454   \textcolor{comment}{! 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   \textcolor{comment}{! 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   \textcolor{keywordflow}{if} ( exp1 > exp2 ) \textcolor{keywordflow}{then}
462     u0\_l = 3.0 * u1 - 2.0 * u0\_r
463 \textcolor{keywordflow}{  endif}
464 
465   \textcolor{keywordflow}{if} ( exp1 < -exp2 ) \textcolor{keywordflow}{then}
466     u0\_r = 3.0 * u1 - 2.0 * u0\_l
467 \textcolor{keywordflow}{  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   \textcolor{comment}{! 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 
\end{DoxyCode}
\mbox{\Hypertarget{namespacepqm__functions_a1e0e86d8470dd334b9cd676f511f6720}\label{namespacepqm__functions_a1e0e86d8470dd334b9cd676f511f6720}} 
\index{pqm\+\_\+functions@{pqm\+\_\+functions}!pqm\+\_\+boundary\+\_\+extrapolation\+\_\+v1@{pqm\+\_\+boundary\+\_\+extrapolation\+\_\+v1}}
\index{pqm\+\_\+boundary\+\_\+extrapolation\+\_\+v1@{pqm\+\_\+boundary\+\_\+extrapolation\+\_\+v1}!pqm\+\_\+functions@{pqm\+\_\+functions}}
\subsubsection{\texorpdfstring{pqm\+\_\+boundary\+\_\+extrapolation\+\_\+v1()}{pqm\_boundary\_extrapolation\_v1()}}
{\footnotesize\ttfamily subroutine, public pqm\+\_\+functions\+::pqm\+\_\+boundary\+\_\+extrapolation\+\_\+v1 (\begin{DoxyParamCaption}\item[{integer, intent(in)}]{N,  }\item[{real, dimension(\+:), intent(in)}]{h,  }\item[{real, dimension(\+:), intent(in)}]{u,  }\item[{real, dimension(\+:,\+:), intent(inout)}]{edge\+\_\+values,  }\item[{real, dimension(\+:,\+:), intent(inout)}]{edge\+\_\+slopes,  }\item[{real, dimension(\+:,\+:), intent(inout)}]{ppoly\+\_\+coef,  }\item[{real, intent(in), optional}]{h\+\_\+neglect }\end{DoxyParamCaption})}



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 P\+PM 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 \textquotesingle{}u\textquotesingle{} is equal to the number of cells defining \textquotesingle{}grid\textquotesingle{} and \textquotesingle{}ppoly\textquotesingle{}. No consistency check is performed here.


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em n} & Number of cells\\
\hline
\mbox{\tt in}  & {\em h} & cell widths (size N) \mbox{[}H\mbox{]}\\
\hline
\mbox{\tt in}  & {\em u} & cell averages (size N) \mbox{[}A\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em edge\+\_\+values} & Edge value of polynomial \mbox{[}A\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em edge\+\_\+slopes} & Edge slope of polynomial \mbox{[}A H-\/1\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em ppoly\+\_\+coef} & Coefficients of polynomial, mainly \mbox{[}A\mbox{]}\\
\hline
\mbox{\tt in}  & {\em h\+\_\+neglect} & A negligibly small width for the purpose of cell reconstructions \mbox{[}H\mbox{]} \\
\hline
\end{DoxyParams}


Definition at line 502 of file P\+Q\+M\+\_\+functions.\+F90.


\begin{DoxyCode}
502   \textcolor{keywordtype}{integer},              \textcolor{keywordtype}{intent(in)}    :: N\textcolor{comment}{ !< Number of cells}
503   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:)},   \textcolor{keywordtype}{intent(in)}    :: h\textcolor{comment}{ !< cell widths (size N) [H]}
504   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:)},   \textcolor{keywordtype}{intent(in)}    :: u\textcolor{comment}{ !< cell averages (size N) [A]}
505   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:,:)}, \textcolor{keywordtype}{intent(inout)} :: edge\_values\textcolor{comment}{    !< Edge value of polynomial [A]}
506   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:,:)}, \textcolor{keywordtype}{intent(inout)} :: edge\_slopes\textcolor{comment}{    !< Edge slope of polynomial [A H-1]}
507   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:,:)}, \textcolor{keywordtype}{intent(inout)} :: ppoly\_coef\textcolor{comment}{ !< Coefficients of polynomial, mainly [A]}
508   \textcolor{keywordtype}{real},       \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}    :: h\_neglect\textcolor{comment}{  !< A negligibly small width for}
509 \textcolor{comment}{                                           !! the purpose of cell reconstructions [H]}
510   \textcolor{comment}{! Local variables}
511   \textcolor{keywordtype}{integer} :: i0, i1
512   \textcolor{keywordtype}{integer} :: inflexion\_l
513   \textcolor{keywordtype}{integer} :: inflexion\_r
514   \textcolor{keywordtype}{real}    :: u0, u1, um
515   \textcolor{keywordtype}{real}    :: h0, h1
516   \textcolor{keywordtype}{real}    :: a, b, c, d, e
517   \textcolor{keywordtype}{real}    :: ar, br, beta
518   \textcolor{keywordtype}{real}    :: u0\_l, u0\_r
519   \textcolor{keywordtype}{real}    :: u1\_l, u1\_r
520   \textcolor{keywordtype}{real}    :: u\_plm
521   \textcolor{keywordtype}{real}    :: slope
522   \textcolor{keywordtype}{real}    :: alpha1, alpha2, alpha3
523   \textcolor{keywordtype}{real}    :: rho, sqrt\_rho
524   \textcolor{keywordtype}{real}    :: gradient1, gradient2
525   \textcolor{keywordtype}{real}    :: x1, x2
526   \textcolor{keywordtype}{real}    :: hNeglect
527 
528   hneglect = hneglect\_dflt ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(h\_neglect)) hneglect = h\_neglect
529 
530   \textcolor{comment}{! ----- 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   \textcolor{comment}{! Compute real slope and express it w.r.t. local coordinate system}
540   \textcolor{comment}{! within boundary cell}
541   slope = 2.0 * ( u1 - u0 ) / ( ( h0 + h1 ) + hneglect )
542   slope = slope * h0
543 
544   \textcolor{comment}{! The right edge value and slope of the boundary cell are taken to be the}
545   \textcolor{comment}{! 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          \textcolor{comment}{! edge value}
550   u1\_r = b / (h1 + hneglect) \textcolor{comment}{! edge slope (w.r.t. global coord.)}
551 
552   \textcolor{comment}{! Compute coefficient for rational function based on mean and right}
553   \textcolor{comment}{! edge value and slope}
554   \textcolor{keywordflow}{if} (u1\_r /= 0.) \textcolor{keywordflow}{then} \textcolor{comment}{! HACK by AJA}
555     beta = 2.0 * ( u0\_r - um ) / ( (h0 + hneglect)*u1\_r) - 1.0
556   \textcolor{keywordflow}{else}
557     beta = 0.
558 \textcolor{keywordflow}{  endif} \textcolor{comment}{! HACK by AJA}
559   br = u0\_r + beta*u0\_r - um
560   ar = um + beta*um - br
561 
562   \textcolor{comment}{! Left edge value estimate based on rational function}
563   u0\_l = ar
564 
565   \textcolor{comment}{! Edge value estimate based on PLM}
566   u\_plm = um - 0.5 * slope
567 
568   \textcolor{comment}{! Check whether the left edge value is bounded by the mean and}
569   \textcolor{comment}{! the PLM edge value. If so, keep it and compute left edge slope}
570   \textcolor{comment}{! based on the rational function. If not, keep the PLM edge value and}
571   \textcolor{comment}{! compute corresponding slope.}
572   \textcolor{keywordflow}{if} ( abs(um-u0\_l) < abs(um-u\_plm) ) \textcolor{keywordflow}{then}
573     u1\_l = 2.0 * ( br - ar*beta)
574     u1\_l = u1\_l / (h0 + hneglect)
575   \textcolor{keywordflow}{else}
576     u0\_l = u\_plm
577     u1\_l = slope / (h0 + hneglect)
578 \textcolor{keywordflow}{  endif}
579 
580   \textcolor{comment}{! 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   \textcolor{comment}{! Check whether inflexion points exist. If so, transform the quartic}
596   \textcolor{comment}{! so that both inflexion points coalesce on the left edge.}
597   \textcolor{keywordflow}{if} (( alpha1 /= 0.0 ) .and. ( rho >= 0.0 )) \textcolor{keywordflow}{then}
598 
599     sqrt\_rho = sqrt( rho )
600 
601     x1 = 0.5 * ( - alpha2 - sqrt\_rho ) / alpha1
602     \textcolor{keywordflow}{if} ( (x1 > 0.0) .and. (x1 < 1.0) ) \textcolor{keywordflow}{then}
603       gradient1 = 4.0 * e * (x1**3) + 3.0 * d * (x1**2) + 2.0 * c * x1 + b
604       \textcolor{keywordflow}{if} ( gradient1 * slope < 0.0 ) \textcolor{keywordflow}{then}
605         inflexion\_l = 1
606 \textcolor{keywordflow}{      endif}
607 \textcolor{keywordflow}{    endif}
608 
609     x2 = 0.5 * ( - alpha2 + sqrt\_rho ) / alpha1
610     \textcolor{keywordflow}{if} ( (x2 > 0.0) .and. (x2 < 1.0) ) \textcolor{keywordflow}{then}
611       gradient2 = 4.0 * e * (x2**3) + 3.0 * d * (x2**2) + 2.0 * c * x2 + b
612       \textcolor{keywordflow}{if} ( gradient2 * slope < 0.0 ) \textcolor{keywordflow}{then}
613         inflexion\_l = 1
614 \textcolor{keywordflow}{      endif}
615 \textcolor{keywordflow}{    endif}
616 
617 \textcolor{keywordflow}{  endif}
618 
619   \textcolor{keywordflow}{if} (( alpha1 == 0.0 ) .and. ( alpha2 /= 0.0 )) \textcolor{keywordflow}{then}
620 
621     x1 = - alpha3 / alpha2
622     \textcolor{keywordflow}{if} ( (x1 >= 0.0) .and. (x1 <= 1.0) ) \textcolor{keywordflow}{then}
623       gradient1 = 3.0 * d * (x1**2) + 2.0 * c * x1 + b
624       \textcolor{keywordflow}{if} ( gradient1 * slope < 0.0 ) \textcolor{keywordflow}{then}
625         inflexion\_l = 1
626 \textcolor{keywordflow}{      endif}
627 \textcolor{keywordflow}{    endif}
628 
629 \textcolor{keywordflow}{  endif}
630 
631   \textcolor{keywordflow}{if} ( inflexion\_l == 1 ) \textcolor{keywordflow}{then}
632 
633     \textcolor{comment}{! We modify the edge slopes so that both inflexion points}
634     \textcolor{comment}{! 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     \textcolor{comment}{! One of the modified slopes might be inconsistent. When that happens,}
639     \textcolor{comment}{! the inconsistent slope is set equal to zero and the opposite edge value}
640     \textcolor{comment}{! and edge slope are modified in compliance with the fact that both}
641     \textcolor{comment}{! inflexion points must still be located on the left edge}
642     \textcolor{keywordflow}{if} ( u1\_l * slope < 0.0 ) \textcolor{keywordflow}{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     \textcolor{keywordflow}{elseif} ( u1\_r * slope < 0.0 ) \textcolor{keywordflow}{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 \textcolor{keywordflow}{    endif}
655 
656 \textcolor{keywordflow}{  endif}
657 
658   \textcolor{comment}{! 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     \textcolor{comment}{! 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   \textcolor{comment}{! ----- 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   \textcolor{comment}{! Compute real slope and express it w.r.t. local coordinate system}
687   \textcolor{comment}{! within boundary cell}
688   slope = 2.0 * ( u1 - u0 ) / ( h0 + h1 )
689   slope = slope * h1
690 
691   \textcolor{comment}{! The left edge value and slope of the boundary cell are taken to be the}
692   \textcolor{comment}{! 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                  \textcolor{comment}{! edge value}
699   u1\_l = (b + 2*c + 3*d + 4*e) / h0         \textcolor{comment}{! edge slope (w.r.t. global coord.)}
700 
701   \textcolor{comment}{! Compute coefficient for rational function based on mean and left}
702   \textcolor{comment}{! edge value and slope}
703   \textcolor{keywordflow}{if} (um-u0\_l /= 0.) \textcolor{keywordflow}{then} \textcolor{comment}{! HACK by AJA}
704     beta = 0.5*h1*u1\_l / (um-u0\_l) - 1.0
705   \textcolor{keywordflow}{else}
706     beta = 0.
707 \textcolor{keywordflow}{  endif} \textcolor{comment}{! HACK by AJA}
708   br = beta*um + um - u0\_l
709   ar = u0\_l
710 
711   \textcolor{comment}{! Right edge value estimate based on rational function}
712   \textcolor{keywordflow}{if} (1+beta /= 0.) \textcolor{keywordflow}{then} \textcolor{comment}{! HACK by AJA}
713     u0\_r = (ar + 2*br + beta*br ) / ((1+beta)*(1+beta))
714   \textcolor{keywordflow}{else}
715     u0\_r = um + 0.5 * slope \textcolor{comment}{! PLM}
716 \textcolor{keywordflow}{  endif} \textcolor{comment}{! HACK by AJA}
717 
718   \textcolor{comment}{! Right edge value estimate based on PLM}
719   u\_plm = um + 0.5 * slope
720 
721   \textcolor{comment}{! Check whether the right edge value is bounded by the mean and}
722   \textcolor{comment}{! the PLM edge value. If so, keep it and compute right edge slope}
723   \textcolor{comment}{! based on the rational function. If not, keep the PLM edge value and}
724   \textcolor{comment}{! compute corresponding slope.}
725   \textcolor{keywordflow}{if} ( abs(um-u0\_r) < abs(um-u\_plm) ) \textcolor{keywordflow}{then}
726     u1\_r = 2.0 * ( br - ar*beta ) / ( (1+beta)*(1+beta)*(1+beta) )
727     u1\_r = u1\_r / h1
728   \textcolor{keywordflow}{else}
729     u0\_r = u\_plm
730     u1\_r = slope / h1
731 \textcolor{keywordflow}{  endif}
732 
733   \textcolor{comment}{! 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   \textcolor{comment}{! Check whether inflexion points exist. If so, transform the quartic}
749   \textcolor{comment}{! so that both inflexion points coalesce on the right edge.}
750   \textcolor{keywordflow}{if} (( alpha1 /= 0.0 ) .and. ( rho >= 0.0 )) \textcolor{keywordflow}{then}
751 
752     sqrt\_rho = sqrt( rho )
753 
754     x1 = 0.5 * ( - alpha2 - sqrt\_rho ) / alpha1
755     \textcolor{keywordflow}{if} ( (x1 > 0.0) .and. (x1 < 1.0) ) \textcolor{keywordflow}{then}
756       gradient1 = 4.0 * e * (x1**3) + 3.0 * d * (x1**2) + 2.0 * c * x1 + b
757       \textcolor{keywordflow}{if} ( gradient1 * slope < 0.0 ) \textcolor{keywordflow}{then}
758         inflexion\_r = 1
759 \textcolor{keywordflow}{      endif}
760 \textcolor{keywordflow}{    endif}
761 
762     x2 = 0.5 * ( - alpha2 + sqrt\_rho ) / alpha1
763     \textcolor{keywordflow}{if} ( (x2 > 0.0) .and. (x2 < 1.0) ) \textcolor{keywordflow}{then}
764       gradient2 = 4.0 * e * (x2**3) + 3.0 * d * (x2**2) + 2.0 * c * x2 + b
765       \textcolor{keywordflow}{if} ( gradient2 * slope < 0.0 ) \textcolor{keywordflow}{then}
766         inflexion\_r = 1
767 \textcolor{keywordflow}{      endif}
768 \textcolor{keywordflow}{    endif}
769 
770 \textcolor{keywordflow}{  endif}
771 
772   \textcolor{keywordflow}{if} (( alpha1 == 0.0 ) .and. ( alpha2 /= 0.0 )) \textcolor{keywordflow}{then}
773 
774     x1 = - alpha3 / alpha2
775     \textcolor{keywordflow}{if} ( (x1 >= 0.0) .and. (x1 <= 1.0) ) \textcolor{keywordflow}{then}
776       gradient1 = 3.0 * d * (x1**2) + 2.0 * c * x1 + b
777       \textcolor{keywordflow}{if} ( gradient1 * slope < 0.0 ) \textcolor{keywordflow}{then}
778         inflexion\_r = 1
779 \textcolor{keywordflow}{      endif}
780 \textcolor{keywordflow}{    endif}
781 
782 \textcolor{keywordflow}{  endif}
783 
784   \textcolor{keywordflow}{if} ( inflexion\_r == 1 ) \textcolor{keywordflow}{then}
785 
786     \textcolor{comment}{! We modify the edge slopes so that both inflexion points}
787     \textcolor{comment}{! 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     \textcolor{comment}{! One of the modified slopes might be inconsistent. When that happens,}
792     \textcolor{comment}{! the inconsistent slope is set equal to zero and the opposite edge value}
793     \textcolor{comment}{! and edge slope are modified in compliance with the fact that both}
794     \textcolor{comment}{! inflexion points must still be located on the right edge}
795     \textcolor{keywordflow}{if} ( u1\_l * slope < 0.0 ) \textcolor{keywordflow}{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     \textcolor{keywordflow}{elseif} ( u1\_r * slope < 0.0 ) \textcolor{keywordflow}{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 \textcolor{keywordflow}{    endif}
808 
809 \textcolor{keywordflow}{  endif}
810 
811   \textcolor{comment}{! 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 
\end{DoxyCode}
\mbox{\Hypertarget{namespacepqm__functions_aab9f3108956943c39ed2d4b675d78021}\label{namespacepqm__functions_aab9f3108956943c39ed2d4b675d78021}} 
\index{pqm\+\_\+functions@{pqm\+\_\+functions}!pqm\+\_\+limiter@{pqm\+\_\+limiter}}
\index{pqm\+\_\+limiter@{pqm\+\_\+limiter}!pqm\+\_\+functions@{pqm\+\_\+functions}}
\subsubsection{\texorpdfstring{pqm\+\_\+limiter()}{pqm\_limiter()}}
{\footnotesize\ttfamily subroutine pqm\+\_\+functions\+::pqm\+\_\+limiter (\begin{DoxyParamCaption}\item[{integer, intent(in)}]{N,  }\item[{real, dimension(\+:), intent(in)}]{h,  }\item[{real, dimension(\+:), intent(in)}]{u,  }\item[{real, dimension(\+:,\+:), intent(inout)}]{edge\+\_\+values,  }\item[{real, dimension(\+:,\+:), intent(inout)}]{edge\+\_\+slopes,  }\item[{real, intent(in), optional}]{h\+\_\+neglect,  }\item[{logical, intent(in), optional}]{answers\+\_\+2018 }\end{DoxyParamCaption})\hspace{0.3cm}{\ttfamily [private]}}



Limit the piecewise quartic method reconstruction. 

Standard P\+QM limiter (White \& Adcroft, J\+CP 2008).

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


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em n} & Number of cells\\
\hline
\mbox{\tt in}  & {\em h} & cell widths (size N) \mbox{[}H\mbox{]}\\
\hline
\mbox{\tt in}  & {\em u} & cell average properties (size N) \mbox{[}A\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em edge\+\_\+values} & Potentially modified edge values \mbox{[}A\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em edge\+\_\+slopes} & Potentially modified edge slopes \mbox{[}A H-\/1\mbox{]}\\
\hline
\mbox{\tt in}  & {\em h\+\_\+neglect} & A negligibly small width for the purpose of cell reconstructions \mbox{[}H\mbox{]}\\
\hline
\mbox{\tt in}  & {\em answers\+\_\+2018} & If true use older, less acccurate expressions. \\
\hline
\end{DoxyParams}


Definition at line 76 of file P\+Q\+M\+\_\+functions.\+F90.


\begin{DoxyCode}
76   \textcolor{keywordtype}{integer},              \textcolor{keywordtype}{intent(in)}    :: N\textcolor{comment}{ !< Number of cells}
77   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:)},   \textcolor{keywordtype}{intent(in)}    :: h\textcolor{comment}{ !< cell widths (size N) [H]}
78   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:)},   \textcolor{keywordtype}{intent(in)}    :: u\textcolor{comment}{ !< cell average properties (size N) [A]}
79   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:,:)}, \textcolor{keywordtype}{intent(inout)} :: edge\_values\textcolor{comment}{ !< Potentially modified edge values [A]}
80   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:,:)}, \textcolor{keywordtype}{intent(inout)} :: edge\_slopes\textcolor{comment}{ !< Potentially modified edge slopes [A H-1]}
81   \textcolor{keywordtype}{real},       \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}    :: h\_neglect\textcolor{comment}{ !< A negligibly small width for}
82 \textcolor{comment}{                                           !! the purpose of cell reconstructions [H]}
83   \textcolor{keywordtype}{logical},    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}    :: answers\_2018\textcolor{comment}{ !< If true use older, less acccurate expressions.}
84 
85   \textcolor{comment}{! Local variables}
86   \textcolor{keywordtype}{integer} :: k            \textcolor{comment}{! loop index}
87   \textcolor{keywordtype}{integer} :: inflexion\_l
88   \textcolor{keywordtype}{integer} :: inflexion\_r
89   \textcolor{keywordtype}{real}    :: u0\_l, u0\_r     \textcolor{comment}{! edge values [A]}
90   \textcolor{keywordtype}{real}    :: u1\_l, u1\_r     \textcolor{comment}{! edge slopes [A H-1]}
91   \textcolor{keywordtype}{real}    :: u\_l, u\_c, u\_r  \textcolor{comment}{! left, center and right cell averages [A]}
92   \textcolor{keywordtype}{real}    :: h\_l, h\_c, h\_r  \textcolor{comment}{! left, center and right cell widths [H]}
93   \textcolor{keywordtype}{real}    :: sigma\_l, sigma\_c, sigma\_r \textcolor{comment}{! left, center and right van Leer slopes}
94   \textcolor{keywordtype}{real}    :: slope          \textcolor{comment}{! retained PLM slope}
95   \textcolor{keywordtype}{real}    :: a, b, c, d, e
96   \textcolor{keywordtype}{real}    :: alpha1, alpha2, alpha3
97   \textcolor{keywordtype}{real}    :: rho, sqrt\_rho
98   \textcolor{keywordtype}{real}    :: gradient1, gradient2
99   \textcolor{keywordtype}{real}    :: x1, x2
100   \textcolor{keywordtype}{real}    :: hNeglect
101 
102   hneglect = hneglect\_dflt ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(h\_neglect)) hneglect = h\_neglect
103 
104   \textcolor{comment}{! Bound edge values}
105   \textcolor{keyword}{call }bound\_edge\_values( n, h, u, edge\_values, hneglect, answers\_2018 )
106 
107   \textcolor{comment}{! Make discontinuous edge values monotonic (thru averaging)}
108   \textcolor{keyword}{call }check\_discontinuous\_edge\_values( n, u, edge\_values )
109 
110   \textcolor{comment}{! Loop on interior cells to apply the PQM limiter}
111   \textcolor{keywordflow}{do} k = 2,n-1
112 
113     \textcolor{comment}{!if ( h(k) < 1.0 ) cycle}
114 
115     inflexion\_l = 0
116     inflexion\_r = 0
117 
118     \textcolor{comment}{! 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     \textcolor{comment}{! Get cell widths and cell averages (boundary cells are assumed to}
125     \textcolor{comment}{! 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     \textcolor{comment}{! 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     \textcolor{keywordflow}{if} ( (sigma\_l * sigma\_r) > 0.0 ) \textcolor{keywordflow}{then}
139       slope = sign( min(abs(sigma\_l),abs(sigma\_c),abs(sigma\_r)), sigma\_c )
140     \textcolor{keywordflow}{else}
141       slope = 0.0
142 \textcolor{keywordflow}{    endif}
143 
144     \textcolor{comment}{! If one of the slopes has the wrong sign compared with the}
145     \textcolor{comment}{! limited PLM slope, it is set equal to the limited PLM slope}
146     \textcolor{keywordflow}{if} ( u1\_l*slope <= 0.0 ) u1\_l = slope
147     \textcolor{keywordflow}{if} ( u1\_r*slope <= 0.0 ) u1\_r = slope
148 
149     \textcolor{comment}{! Local extremum --> flatten}
150     \textcolor{keywordflow}{if} ( (u0\_r - u\_c) * (u\_c - u0\_l) <= 0.0) \textcolor{keywordflow}{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 \textcolor{keywordflow}{    endif}
158 
159     \textcolor{comment}{! Edge values are bounded and averaged when discontinuous and not}
160     \textcolor{comment}{! monotonic, edge slopes are consistent and the cell is not an extremum.}
161     \textcolor{comment}{! We now need to check and encorce the monotonicity of the quartic within}
162     \textcolor{comment}{! the cell}
163     \textcolor{keywordflow}{if} ( (inflexion\_l == 0) .AND. (inflexion\_r == 0) ) \textcolor{keywordflow}{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       \textcolor{comment}{! Determine the coefficients of the second derivative}
172       \textcolor{comment}{! 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       \textcolor{comment}{! Check whether inflexion points exist}
180       \textcolor{keywordflow}{if} (( alpha1 /= 0.0 ) .and. ( rho >= 0.0 )) \textcolor{keywordflow}{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         \textcolor{comment}{! Check whether both inflexion points lie in [0,1]}
188         \textcolor{keywordflow}{if} ( (x1 >= 0.0) .AND. (x1 <= 1.0) .AND. &
189              (x2 >= 0.0) .AND. (x2 <= 1.0) ) \textcolor{keywordflow}{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           \textcolor{comment}{! Check whether one of the gradients is inconsistent}
195           \textcolor{keywordflow}{if} ( (gradient1 * slope < 0.0) .OR. &
196                (gradient2 * slope < 0.0) ) \textcolor{keywordflow}{then}
197             \textcolor{comment}{! Decide where to collapse inflexion points}
198             \textcolor{comment}{! (depends on one-sided slopes)}
199             \textcolor{keywordflow}{if} ( abs(sigma\_l) < abs(sigma\_r) ) \textcolor{keywordflow}{then}
200               inflexion\_l = 1
201             \textcolor{keywordflow}{else}
202               inflexion\_r = 1
203 \textcolor{keywordflow}{            endif}
204 \textcolor{keywordflow}{          endif}
205 
206         \textcolor{comment}{! If both x1 and x2 do not lie in [0,1], check whether}
207         \textcolor{comment}{! only x1 lies in [0,1]}
208         \textcolor{keywordflow}{elseif} ( (x1 >= 0.0) .AND. (x1 <= 1.0) ) \textcolor{keywordflow}{then}
209 
210           gradient1 = 4.0 * e * (x1**3) + 3.0 * d * (x1**2) + 2.0 * c * x1 + b
211 
212           \textcolor{comment}{! Check whether the gradient is inconsistent}
213           \textcolor{keywordflow}{if} ( gradient1 * slope < 0.0 ) \textcolor{keywordflow}{then}
214             \textcolor{comment}{! Decide where to collapse inflexion points}
215             \textcolor{comment}{! (depends on one-sided slopes)}
216             \textcolor{keywordflow}{if} ( abs(sigma\_l) < abs(sigma\_r) ) \textcolor{keywordflow}{then}
217               inflexion\_l = 1
218             \textcolor{keywordflow}{else}
219               inflexion\_r = 1
220 \textcolor{keywordflow}{            endif}
221 \textcolor{keywordflow}{          endif}
222 
223         \textcolor{comment}{! If x1 does not lie in [0,1], check whether x2 lies in [0,1]}
224         \textcolor{keywordflow}{elseif} ( (x2 >= 0.0) .AND. (x2 <= 1.0) ) \textcolor{keywordflow}{then}
225 
226           gradient2 = 4.0 * e * (x2**3) + 3.0 * d * (x2**2) + 2.0 * c * x2 + b
227 
228           \textcolor{comment}{! Check whether the gradient is inconsistent}
229           \textcolor{keywordflow}{if} ( gradient2 * slope < 0.0 ) \textcolor{keywordflow}{then}
230             \textcolor{comment}{! Decide where to collapse inflexion points}
231             \textcolor{comment}{! (depends on one-sided slopes)}
232             \textcolor{keywordflow}{if} ( abs(sigma\_l) < abs(sigma\_r) ) \textcolor{keywordflow}{then}
233               inflexion\_l = 1
234             \textcolor{keywordflow}{else}
235               inflexion\_r = 1
236 \textcolor{keywordflow}{            endif}
237 \textcolor{keywordflow}{          endif}
238 
239 \textcolor{keywordflow}{        endif} \textcolor{comment}{! end checking where the inflexion points lie}
240 
241 \textcolor{keywordflow}{      endif} \textcolor{comment}{! end checking if alpha1 != 0 AND rho >= 0}
242 
243       \textcolor{comment}{! If alpha1 is zero, the second derivative of the quartic reduces}
244       \textcolor{comment}{! to a straight line}
245       \textcolor{keywordflow}{if} (( alpha1 == 0.0 ) .and. ( alpha2 /= 0.0 )) \textcolor{keywordflow}{then}
246 
247           x1 = - alpha3 / alpha2
248           \textcolor{keywordflow}{if} ( (x1 >= 0.0) .AND. (x1 <= 1.0) ) \textcolor{keywordflow}{then}
249 
250             gradient1 = 4.0 * e * (x1**3) + 3.0 * d * (x1**2) + 2.0 * c * x1 + b
251 
252             \textcolor{comment}{! Check whether the gradient is inconsistent}
253             \textcolor{keywordflow}{if} ( gradient1 * slope < 0.0 ) \textcolor{keywordflow}{then}
254               \textcolor{comment}{! Decide where to collapse inflexion points}
255               \textcolor{comment}{! (depends on one-sided slopes)}
256               \textcolor{keywordflow}{if} ( abs(sigma\_l) < abs(sigma\_r) ) \textcolor{keywordflow}{then}
257                 inflexion\_l = 1
258               \textcolor{keywordflow}{else}
259                 inflexion\_r = 1
260 \textcolor{keywordflow}{              endif}
261 \textcolor{keywordflow}{            endif} \textcolor{comment}{! check slope consistency}
262 
263 \textcolor{keywordflow}{          endif}
264 
265 \textcolor{keywordflow}{      endif} \textcolor{comment}{! end check whether we can find the root of the straight line}
266 
267 \textcolor{keywordflow}{    endif} \textcolor{comment}{! end checking whether to shift inflexion points}
268 
269     \textcolor{comment}{! At this point, we know onto which edge to shift inflexion points}
270     \textcolor{keywordflow}{if} ( inflexion\_l == 1 ) \textcolor{keywordflow}{then}
271 
272       \textcolor{comment}{! We modify the edge slopes so that both inflexion points}
273       \textcolor{comment}{! 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       \textcolor{comment}{! One of the modified slopes might be inconsistent. When that happens,}
278       \textcolor{comment}{! the inconsistent slope is set equal to zero and the opposite edge value}
279       \textcolor{comment}{! and edge slope are modified in compliance with the fact that both}
280       \textcolor{comment}{! inflexion points must still be located on the left edge}
281       \textcolor{keywordflow}{if} ( u1\_l * slope < 0.0 ) \textcolor{keywordflow}{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       \textcolor{keywordflow}{elseif} ( u1\_r * slope < 0.0 ) \textcolor{keywordflow}{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 \textcolor{keywordflow}{      endif}
294 
295     \textcolor{keywordflow}{elseif} ( inflexion\_r == 1 ) \textcolor{keywordflow}{then}
296 
297       \textcolor{comment}{! We modify the edge slopes so that both inflexion points}
298       \textcolor{comment}{! 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       \textcolor{comment}{! One of the modified slopes might be inconsistent. When that happens,}
303       \textcolor{comment}{! the inconsistent slope is set equal to zero and the opposite edge value}
304       \textcolor{comment}{! and edge slope are modified in compliance with the fact that both}
305       \textcolor{comment}{! inflexion points must still be located on the right edge}
306       \textcolor{keywordflow}{if} ( u1\_l * slope < 0.0 ) \textcolor{keywordflow}{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       \textcolor{keywordflow}{elseif} ( u1\_r * slope < 0.0 ) \textcolor{keywordflow}{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 \textcolor{keywordflow}{      endif}
319 
320 \textcolor{keywordflow}{    endif} \textcolor{comment}{! clause to check where to collapse inflexion points}
321 
322     \textcolor{comment}{! 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 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! end loop on interior cells}
329 
330   \textcolor{comment}{! 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 
\end{DoxyCode}
\mbox{\Hypertarget{namespacepqm__functions_afa6f7b5430011f03c428b329c5f42fae}\label{namespacepqm__functions_afa6f7b5430011f03c428b329c5f42fae}} 
\index{pqm\+\_\+functions@{pqm\+\_\+functions}!pqm\+\_\+reconstruction@{pqm\+\_\+reconstruction}}
\index{pqm\+\_\+reconstruction@{pqm\+\_\+reconstruction}!pqm\+\_\+functions@{pqm\+\_\+functions}}
\subsubsection{\texorpdfstring{pqm\+\_\+reconstruction()}{pqm\_reconstruction()}}
{\footnotesize\ttfamily subroutine, public pqm\+\_\+functions\+::pqm\+\_\+reconstruction (\begin{DoxyParamCaption}\item[{integer, intent(in)}]{N,  }\item[{real, dimension(\+:), intent(in)}]{h,  }\item[{real, dimension(\+:), intent(in)}]{u,  }\item[{real, dimension(\+:,\+:), intent(inout)}]{edge\+\_\+values,  }\item[{real, dimension(\+:,\+:), intent(inout)}]{edge\+\_\+slopes,  }\item[{real, dimension(\+:,\+:), intent(inout)}]{ppoly\+\_\+coef,  }\item[{real, intent(in), optional}]{h\+\_\+neglect,  }\item[{logical, intent(in), optional}]{answers\+\_\+2018 }\end{DoxyParamCaption})}



Reconstruction by quartic polynomials within each cell. 

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


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em n} & Number of cells\\
\hline
\mbox{\tt in}  & {\em h} & cell widths (size N) \mbox{[}H\mbox{]}\\
\hline
\mbox{\tt in}  & {\em u} & cell averages (size N) \mbox{[}A\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em edge\+\_\+values} & Edge value of polynomial \mbox{[}A\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em edge\+\_\+slopes} & Edge slope of polynomial \mbox{[}A H-\/1\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em ppoly\+\_\+coef} & Coefficients of polynomial, mainly \mbox{[}A\mbox{]}\\
\hline
\mbox{\tt in}  & {\em h\+\_\+neglect} & A negligibly small width for the purpose of cell reconstructions \mbox{[}H\mbox{]}\\
\hline
\mbox{\tt in}  & {\em answers\+\_\+2018} & If true use older, less acccurate expressions. \\
\hline
\end{DoxyParams}


Definition at line 21 of file P\+Q\+M\+\_\+functions.\+F90.


\begin{DoxyCode}
21   \textcolor{keywordtype}{integer},              \textcolor{keywordtype}{intent(in)}    :: N\textcolor{comment}{ !< Number of cells}
22   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:)},   \textcolor{keywordtype}{intent(in)}    :: h\textcolor{comment}{ !< cell widths (size N) [H]}
23   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:)},   \textcolor{keywordtype}{intent(in)}    :: u\textcolor{comment}{ !< cell averages (size N) [A]}
24   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:,:)}, \textcolor{keywordtype}{intent(inout)} :: edge\_values\textcolor{comment}{    !< Edge value of polynomial [A]}
25   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:,:)}, \textcolor{keywordtype}{intent(inout)} :: edge\_slopes\textcolor{comment}{    !< Edge slope of polynomial [A H-1]}
26   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(:,:)}, \textcolor{keywordtype}{intent(inout)} :: ppoly\_coef\textcolor{comment}{ !< Coefficients of polynomial, mainly [A]}
27   \textcolor{keywordtype}{real},       \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}    :: h\_neglect\textcolor{comment}{  !< A negligibly small width for}
28 \textcolor{comment}{                                           !! the purpose of cell reconstructions [H]}
29   \textcolor{keywordtype}{logical},    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}    :: answers\_2018\textcolor{comment}{ !< If true use older, less acccurate expressions.}
30 
31   \textcolor{comment}{! Local variables}
32   \textcolor{keywordtype}{integer}   :: k                \textcolor{comment}{! loop index}
33   \textcolor{keywordtype}{real}      :: h\_c              \textcolor{comment}{! cell width}
34   \textcolor{keywordtype}{real}      :: u0\_l, u0\_r       \textcolor{comment}{! edge values (left and right) [A]}
35   \textcolor{keywordtype}{real}      :: u1\_l, u1\_r       \textcolor{comment}{! edge slopes (left and right) [A H-1]}
36   \textcolor{keywordtype}{real}      :: a, b, c, d, e    \textcolor{comment}{! parabola coefficients}
37 
38   \textcolor{comment}{! PQM limiter}
39   \textcolor{keyword}{call }pqm\_limiter( n, h, u, edge\_values, edge\_slopes, h\_neglect, answers\_2018 )
40 
41   \textcolor{comment}{! Loop on cells to construct the cubic within each cell}
42   \textcolor{keywordflow}{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     \textcolor{comment}{! 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 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! end loop on cells}
66 
\end{DoxyCode}
