From 5bde4b25c49461e9ead4f1f5bc1eaff449ad6ce5 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Fri, 9 Aug 2024 12:26:32 -0400 Subject: [PATCH] Force bounds on x argument --- src/ALE/Recon1d_PLM_CW.F90 | 6 ++++-- src/ALE/Recon1d_PPM_H4_2019.F90 | 12 ++++++++---- 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/src/ALE/Recon1d_PLM_CW.F90 b/src/ALE/Recon1d_PLM_CW.F90 index 7b70ce0006..760b0bc902 100644 --- a/src/ALE/Recon1d_PLM_CW.F90 +++ b/src/ALE/Recon1d_PLM_CW.F90 @@ -168,15 +168,17 @@ real function f(this, k, x) class(PLM_CW), intent(in) :: this !< This reconstruction integer, intent(in) :: k !< Cell number real, intent(in) :: x !< Non-dimensional position within element [nondim] + real :: xc ! Bounded version of x [nondim] real :: du ! Difference across cell [A] real :: u_a, u_b ! Two estimate of f [A] du = this%ur(k) - this%ul(k) + xc = max( 0., min( 1., x ) ) ! This expression for u_a can overshoot u_r but is good for x<<1 - u_a = this%ul(k) + du * x + u_a = this%ul(k) + du * xc ! This expression for u_b can overshoot u_l but is good for 1-x<<1 - u_b = this%ur(k) + du * ( x - 1. ) + u_b = this%ur(k) + du * ( xc - 1. ) ! Since u_a and u_b are both bounded, this will perserve uniformity f = 0.5 * ( u_a + u_b ) diff --git a/src/ALE/Recon1d_PPM_H4_2019.F90 b/src/ALE/Recon1d_PPM_H4_2019.F90 index e031fd5d01..7314ad4acb 100644 --- a/src/ALE/Recon1d_PPM_H4_2019.F90 +++ b/src/ALE/Recon1d_PPM_H4_2019.F90 @@ -119,6 +119,7 @@ real function f(this, k, x) class(PPM_H4_2019), intent(in) :: this !< This reconstruction integer, intent(in) :: k !< Cell number real, intent(in) :: x !< Non-dimensional position within element [nondim] + real :: xc ! Bounded version of x [nondim] real :: du ! Difference across cell [A] real :: a6 ! Collela and Woordward curvature parameter [A] real :: u_a, u_b ! Two estimate of f [A] @@ -127,15 +128,16 @@ real function f(this, k, x) du = this%ur(k) - this%ul(k) a6 = 3.0 * ( ( this%u_mean(k) - this%ul(k) ) + ( this%u_mean(k) - this%ur(k) ) ) - lmx = 1.0 - x + xc = max( 0., min( 1., x ) ) + lmx = 1.0 - xc ! This expression for u_a can overshoot u_r but is good for x<<1 - u_a = this%ul(k) + x * ( du + a6 * x ) + u_a = this%ul(k) + xc * ( du + a6 * xc ) ! This expression for u_b can overshoot u_l but is good for 1-x<<1 u_b = this%ur(k) + lmx * ( - du + a6 * lmx ) ! Since u_a and u_b are both side-bounded, using weights=0 or 1 will preserve uniformity - wb = 0.5 + sign(0.5, x - 0.5 ) ! = 1 @ x=0, = 0 @ x=1 + wb = 0.5 + sign(0.5, xc - 0.5 ) ! = 1 @ x=0, = 0 @ x=1 f = ( ( 1. - wb ) * u_a ) + ( wb * u_b ) end function f @@ -145,13 +147,15 @@ real function dfdx(this, k, x) class(PPM_H4_2019), intent(in) :: this !< This reconstruction integer, intent(in) :: k !< Cell number real, intent(in) :: x !< Non-dimensional position within element [nondim] + real :: xc ! Bounded version of x [nondim] real :: du ! Difference across cell [A] real :: a6 ! Collela and Woordward curvature parameter [A] du = this%ur(k) - this%ul(k) a6 = 3.0 * ( ( this%u_mean(k) - this%ul(k) ) + ( this%u_mean(k) - this%ur(k) ) ) + xc = max( 0., min( 1., x ) ) - dfdx = du + a6 * ( 2.0 * x - 1.0 ) + dfdx = du + a6 * ( 2.0 * xc - 1.0 ) end function dfdx