\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}