Skip to content

Commit

Permalink
Force bounds on x argument
Browse files Browse the repository at this point in the history
  • Loading branch information
adcroft committed Aug 9, 2024
1 parent f260a46 commit 5bde4b2
Show file tree
Hide file tree
Showing 2 changed files with 12 additions and 6 deletions.
6 changes: 4 additions & 2 deletions src/ALE/Recon1d_PLM_CW.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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 )
Expand Down
12 changes: 8 additions & 4 deletions src/ALE/Recon1d_PPM_H4_2019.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand All @@ -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
Expand All @@ -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

Expand Down

0 comments on commit 5bde4b2

Please sign in to comment.