Skip to content

Commit

Permalink
Inlined fns for PPM_H4_2019
Browse files Browse the repository at this point in the history
  • Loading branch information
@adcroft
adcroft committed Oct 4, 2024
1 parent 79ae180 commit 1793985
Showing 1 changed file with 262 additions and 10 deletions.
272 changes: 262 additions & 10 deletions src/ALE/Recon1d_PPM_H4_2019.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,14 +3,12 @@
!! This implementation of PPM ostensbily follows Colella and Woodward, 1984, with cells resorting to PCM for
!! extrema including first and last cells in column. ++++++++++++++++++++++
!! The first and last cells are always limited to PCM.
!! This uses numerical expressions refactored at the beginning of 2019.
module Recon1d_PPM_H4_2019

! This file is part of MOM6. See LICENSE.md for the license.

use Recon1d_type, only : Recon1d, testing
use regrid_edge_values, only : edge_values_explicit_h4
use regrid_edge_values, only : bound_edge_values, check_discontinuous_edge_values
use PPM_functions, only : PPM_reconstruction, PPM_boundary_extrapolation, PPM_monotonicity

implicit none ; private

Expand Down Expand Up @@ -95,15 +93,171 @@ subroutine reconstruct(this, h, u)
real :: u_l, u_r, u_c ! Left, right, and center values [A]
real :: h_l, h_c, h_r ! Thickness of left, center and right cells [H]
real :: h_c0 ! Thickness of center with h_neglect added [H]
integer :: k, n
real :: h0, h1, h2, h3 ! temporary thicknesses [H]
real :: h_min ! A minimal cell width [H]
real :: f1 ! An auxiliary variable [H]
real :: f2 ! An auxiliary variable [A H]
real :: f3 ! An auxiliary variable [H-1]
real :: et1, et2, et3 ! terms the expression for edge values [A H]
real :: I_h12 ! The inverse of the sum of the two central thicknesses [H-1]
real :: I_h012, I_h123 ! Inverses of sums of three successive thicknesses [H-1]
real :: I_den_et2, I_den_et3 ! Inverses of denominators in edge value terms [H-2]
real :: dx ! Difference of successive values of x [H]
real :: f ! value of polynomial at x in arbitrary units [A]
real :: edge_l, edge_r ! Edge values (left and right) [A]
real :: expr1, expr2 ! Temporary expressions [A2]
real :: slope_x_h ! retained PLM slope times half grid step [A]
real :: u0_avg ! avg value at given edge [A]
real, parameter :: hMinFrac = 1.e-5 !< A minimum fraction for min(h)/sum(h) [nondim]
real :: edge_values(this%n,2) ! Edge values [A]
real :: ppoly_coef(this%n,3) ! Polynomial coefficients [A]
real :: ppoly_coef(this%n,3) ! Polynomial coefficients [A]
real :: dz(4) ! A temporary array of limited layer thicknesses [H]
real :: u_tmp(4) ! A temporary array of cell average properties [A]
real :: A(4,4) ! Differences in successive positions raised to various powers,
! in units that vary with the second (j) index as [H^j]
real :: B(4) ! The right hand side of the system to solve for C [A H]
real :: C(4) ! The coefficients of a fit polynomial in units that vary
! with the index (j) as [A H^(j-1)]
integer :: k, n, km1, kp1

n = this%n
call edge_values_explicit_h4( n, h, u, edge_values, h_neglect=this%h_neglect, &
answer_date=20190101 )
call PPM_reconstruction( N, h, u, edge_values, ppoly_coef, this%h_neglect, &
answer_date=20190101 )

! Loop on interior cells
do K = 3, n-1

h0 = h(k-2)
h1 = h(k-1)
h2 = h(k)
h3 = h(k+1)

! Avoid singularities when consecutive pairs of h vanish
if (h0+h1==0.0 .or. h1+h2==0.0 .or. h2+h3==0.0) then
h_min = hMinFrac*max( this%h_neglect, (h0+h1)+(h2+h3) )
h0 = max( h_min, h0 )
h1 = max( h_min, h1 )
h2 = max( h_min, h2 )
h3 = max( h_min, h3 )
endif

I_h12 = 1.0 / (h1+h2)
I_den_et2 = 1.0 / ( ((h0+h1)+h2)*(h0+h1) ) ; I_h012 = (h0+h1) * I_den_et2
I_den_et3 = 1.0 / ( (h1+(h2+h3))*(h2+h3) ) ; I_h123 = (h2+h3) * I_den_et3

et1 = ( 1.0 + (h1 * I_h012 + (h0+h1) * I_h123) ) * I_h12 * (h2*(h2+h3)) * u(k-1) + &
( 1.0 + (h2 * I_h123 + (h2+h3) * I_h012) ) * I_h12 * (h1*(h0+h1)) * u(k)
et2 = ( h1 * (h2*(h2+h3)) * I_den_et2 ) * (u(k-1)-u(k-2))
et3 = ( h2 * (h1*(h0+h1)) * I_den_et3 ) * (u(k) - u(k+1))
edge_values(k,1) = (et1 + (et2 + et3)) / ((h0 + h1) + (h2 + h3))
edge_values(k-1,2) = edge_values(k,1)

enddo ! end loop on interior cells

! Determine first two edge values
do k=1,4 ; dz(k) = max(this%h_neglect, h(k) ) ; u_tmp(k) = u(k) ; enddo
call end_value_h4(dz, u_tmp, C)

! Set the edge values of the first cell
edge_values(1,1) = C(1)
edge_values(1,2) = C(1) + dz(1) * ( C(2) + dz(1) * ( C(3) + dz(1) * C(4) ) )
edge_values(2,1) = edge_values(1,2)

! Determine two edge values of the last cell
do k=1,4 ; dz(k) = max(this%h_neglect, h(n+1-k) ) ; u_tmp(k) = u(n+1-k) ; enddo
call end_value_h4(dz, u_tmp, C)

! Set the last and second to last edge values
edge_values(n,2) = C(1)
edge_values(n,1) = C(1) + dz(1) * ( C(2) + dz(1) * ( C(3) + dz(1) * C(4) ) )
edge_values(n-1,2) = edge_values(n,1)

! Loop on cells to bound edge value
do k = 1, n

! For the sake of bounding boundary edge values, the left neighbor of the left boundary cell
! is assumed to be the same as the left boundary cell and the right neighbor of the right
! boundary cell is assumed to be the same as the right boundary cell. This effectively makes
! boundary cells look like extrema.
km1 = max(1,k-1) ; kp1 = min(k+1,N)

slope_x_h = 0.0
sigma_l = ( u(k) - u(km1) )
sigma_c = ( u(kp1) - u(km1) ) * ( h(k) / ((h(km1) + h(kp1)) + 2.0*h(k)) )
sigma_r = ( u(kp1) - u(k) )

! The limiter is used in the local coordinate system to each cell, so for convenience store
! the slope times a half grid spacing. (See White and Adcroft JCP 2008 Eqs 19 and 20)
if ( (sigma_l * sigma_r) > 0.0 ) &
slope_x_h = sign( min(abs(sigma_l),abs(sigma_c),abs(sigma_r)), sigma_c )

! Limit the edge values
if ( (u(km1)-edge_values(k,1)) * (edge_values(k,1)-u(k)) < 0.0 ) then
edge_values(k,1) = u(k) - sign( min( abs(slope_x_h), abs(edge_values(k,1)-u(k)) ), slope_x_h )
endif

if ( (u(kp1)-edge_values(k,2)) * (edge_values(k,2)-u(k)) < 0.0 ) then
edge_values(k,2) = u(k) + sign( min( abs(slope_x_h), abs(edge_values(k,2)-u(k)) ), slope_x_h )
endif

! Finally bound by neighboring cell means in case of roundoff
edge_values(k,1) = max( min( edge_values(k,1), max(u(km1), u(k)) ), min(u(km1), u(k)) )
edge_values(k,2) = max( min( edge_values(k,2), max(u(kp1), u(k)) ), min(u(kp1), u(k)) )

enddo ! loop on interior edges

do k = 1, n-1
if ( (edge_values(k+1,1) - edge_values(k,2)) * (u(k+1) - u(k)) < 0.0 ) then
u0_avg = 0.5 * ( edge_values(k,2) + edge_values(k+1,1) )
u0_avg = max( min( u0_avg, max(u(k), u(k+1)) ), min(u(k), u(k+1)) )
edge_values(k,2) = u0_avg
edge_values(k+1,1) = u0_avg
endif
enddo ! end loop on interior edges

! Loop on interior cells to apply the standard
! PPM limiter (Colella & Woodward, JCP 84)
do k = 2,N-1

! Get cell averages
u_l = u(k-1)
u_c = u(k)
u_r = u(k+1)

edge_l = edge_values(k,1)
edge_r = edge_values(k,2)

if ( (u_r - u_c)*(u_c - u_l) <= 0.0) then
! Flatten extremum
edge_l = u_c
edge_r = u_c
else
expr1 = 3.0 * (edge_r - edge_l) * ( (u_c - edge_l) + (u_c - edge_r))
expr2 = (edge_r - edge_l) * (edge_r - edge_l)
if ( expr1 > expr2 ) then
! Place extremum at right edge of cell by adjusting left edge value
edge_l = u_c + 2.0 * ( u_c - edge_r )
edge_l = max( min( edge_l, max(u_l, u_c) ), min(u_l, u_c) ) ! In case of round off
elseif ( expr1 < -expr2 ) then
! Place extremum at left edge of cell by adjusting right edge value
edge_r = u_c + 2.0 * ( u_c - edge_l )
edge_r = max( min( edge_r, max(u_r, u_c) ), min(u_r, u_c) ) ! In case of round off
endif
endif
! This checks that the difference in edge values is representable
! and avoids overshoot problems due to round off.
!### The 1.e-60 needs to have units of [A], so this dimensionally inconsistent.
if ( abs( edge_r - edge_l )<max(1.e-60,epsilon(u_c)*abs(u_c)) ) then
edge_l = u_c
edge_r = u_c
endif

edge_values(k,1) = edge_l
edge_values(k,2) = edge_r

enddo ! end loop on interior cells

! PCM within boundary cells
edge_values(1,:) = u(1)
edge_values(N,:) = u(N)

! Store reconstruction
do k = 1, n
Expand All @@ -114,6 +268,104 @@ subroutine reconstruct(this, h, u)

end subroutine reconstruct

!> Determine a one-sided 4th order polynomial fit of u to the data points for the purposes of specifying
!! edge values, as described in the appendix of White and Adcroft JCP 2008.
subroutine end_value_h4(dz, u, Csys)
real, intent(in) :: dz(4) !< The thicknesses of 4 layers, starting at the edge [H].
!! The values of dz must be positive.
real, intent(in) :: u(4) !< The average properties of 4 layers, starting at the edge [A]
real, intent(out) :: Csys(4) !< The four coefficients of a 4th order polynomial fit
!! of u as a function of z [A H-(n-1)]

! Local variables
real :: Wt(3,4) ! The weights of successive u differences in the 4 closed form expressions.
! The units of Wt vary with the second index as [H-(n-1)].
real :: h1, h2, h3, h4 ! Copies of the layer thicknesses [H]
real :: h12, h23, h34 ! Sums of two successive thicknesses [H]
real :: h123, h234 ! Sums of three successive thicknesses [H]
real :: h1234 ! Sums of all four thicknesses [H]
! real :: I_h1 ! The inverse of the a thickness [H-1]
real :: I_h12, I_h23, I_h34 ! The inverses of sums of two thicknesses [H-1]
real :: I_h123, I_h234 ! The inverse of the sum of three thicknesses [H-1]
real :: I_h1234 ! The inverse of the sum of all four thicknesses [H-1]
real :: I_denom ! The inverse of the denominator some expressions [H-3]
real :: I_denB3 ! The inverse of the product of three sums of thicknesses [H-3]
real :: min_frac = 1.0e-6 ! The square of min_frac should be much larger than roundoff [nondim]
real, parameter :: C1_3 = 1.0 / 3.0 ! A rational parameter [nondim]

! if ((dz(1) == dz(2)) .and. (dz(1) == dz(3)) .and. (dz(1) == dz(4))) then
! ! There are simple closed-form expressions in this case
! I_h1 = 0.0 ; if (dz(1) > 0.0) I_h1 = 1.0 / dz(1)
! Csys(1) = u(1) + (-13.0 * (u(2)-u(1)) + 10.0 * (u(3)-u(2)) - 3.0 * (u(4)-u(3))) * (0.25*C1_3)
! Csys(2) = (35.0 * (u(2)-u(1)) - 34.0 * (u(3)-u(2)) + 11.0 * (u(4)-u(3))) * (0.25*C1_3 * I_h1)
! Csys(3) = (-5.0 * (u(2)-u(1)) + 8.0 * (u(3)-u(2)) - 3.0 * (u(4)-u(3))) * (0.25 * I_h1**2)
! Csys(4) = ((u(2)-u(1)) - 2.0 * (u(3)-u(2)) + (u(4)-u(3))) * (0.5*C1_3)
! else

! Express the coefficients as sums of the differences between properties of successive layers.

h1 = dz(1) ; h2 = dz(2) ; h3 = dz(3) ; h4 = dz(4)
! Some of the weights used below are proportional to (h1/(h2+h3))**2 or (h1/(h2+h3))*(h2/(h3+h4))
! so h2 and h3 should be adjusted to ensure that these ratios are not so large that property
! differences at the level of roundoff are amplified to be of order 1.
if ((h2+h3) < min_frac*h1) h3 = min_frac*h1 - h2
if ((h3+h4) < min_frac*h1) h4 = min_frac*h1 - h3

h12 = h1+h2 ; h23 = h2+h3 ; h34 = h3+h4
h123 = h12 + h3 ; h234 = h2 + h34 ; h1234 = h12 + h34
! Find 3 reciprocals with a single division for efficiency.
I_denB3 = 1.0 / (h123 * h12 * h23)
I_h12 = (h123 * h23) * I_denB3
I_h23 = (h12 * h123) * I_denB3
I_h123 = (h12 * h23) * I_denB3
I_denom = 1.0 / ( h1234 * (h234 * h34) )
I_h34 = (h1234 * h234) * I_denom
I_h234 = (h1234 * h34) * I_denom
I_h1234 = (h234 * h34) * I_denom

! Calculation coefficients in the four equations

! The expressions for Csys(3) and Csys(4) come from reducing the 4x4 matrix problem into the following 2x2
! matrix problem, then manipulating the analytic solution to avoid any subtraction and simplifying.
! (C1_3 * h123 * h23) * Csys(3) + (0.25 * h123 * h23 * (h3 + 2.0*h2 + 3.0*h1)) * Csys(4) =
! (u(3)-u(1)) - (u(2)-u(1)) * (h12 + h23) * I_h12
! (C1_3 * ((h23 + h34) * h1234 + h23 * h3)) * Csys(3) +
! (0.25 * ((h1234 + h123 + h12 + h1) * h23 * h3 + (h1234 + h12 + h1) * (h23 + h34) * h1234)) * Csys(4) =
! (u(4)-u(1)) - (u(2)-u(1)) * (h123 + h234) * I_h12
! The final expressions for Csys(1) and Csys(2) were derived by algebraically manipulating the following expressions:
! Csys(1) = (C1_3 * h1 * h12 * Csys(3) + 0.25 * h1 * h12 * (2.0*h1+h2) * Csys(4)) - &
! (h1*I_h12)*(u(2)-u(1)) + u(1)
! Csys(2) = (-2.0*C1_3 * (2.0*h1+h2) * Csys(3) - 0.5 * (h1**2 + h12 * (2.0*h1+h2)) * Csys(4)) + &
! 2.0*I_h12 * (u(2)-u(1))
! These expressions are typically evaluated at x=0 and x=h1, so it is important that these are well behaved
! for these values, suggesting that h1/h23 and h1/h34 should not be allowed to be too large.

Wt(1,1) = -h1 * (I_h1234 + I_h123 + I_h12) ! > -3
Wt(2,1) = h1 * h12 * ( I_h234 * I_h1234 + I_h23 * (I_h234 + I_h123) ) ! < (h1/h234) + (h1/h23)*(2+(h1/h234))
Wt(3,1) = -h1 * h12 * h123 * I_denom ! > -(h1/h34)*(1+(h1/h234))

Wt(1,2) = 2.0 * (I_h12*(1.0 + (h1+h12) * (I_h1234 + I_h123)) + h1 * I_h1234*I_h123) ! < 10/h12
Wt(2,2) = -2.0 * ((h1 * h12 * I_h1234) * (I_h23 * (I_h234 + I_h123)) + & ! > -(10+6*(h1/h234))/h23
(h1+h12) * ( I_h1234*I_h234 + I_h23 * (I_h234 + I_h123) ) )
Wt(3,2) = 2.0 * ((h1+h12) * h123 + h1*h12 ) * I_denom ! < (2+(6*h1/h234)) / h34

Wt(1,3) = -3.0 * I_h12 * I_h123* ( 1.0 + I_h1234 * ((h1+h12)+h123) ) ! > -12 / (h12*h123)
Wt(2,3) = 3.0 * I_h23 * ( I_h123 + I_h1234 * ((h1+h12)+h123) * (I_h123 + I_h234) ) ! < 12 / (h23^2)
Wt(3,3) = -3.0 * ((h1+h12)+h123) * I_denom ! > -9 / (h234*h23)

Wt(1,4) = 4.0 * I_h1234 * I_h123 * I_h12 ! Wt*h1^3 < 4
Wt(2,4) = -4.0 * I_h1234 * (I_h23 * (I_h123 + I_h234)) ! Wt*h1^3 > -4* (h1/h23)*(1+h1/h234)
Wt(3,4) = 4.0 * I_denom ! = 4.0*I_h1234 * I_h234 * I_h34 ! Wt*h1^3 < 4 * (h1/h234)*(h1/h34)

Csys(1) = ((u(1) + Wt(1,1) * (u(2)-u(1))) + Wt(2,1) * (u(3)-u(2))) + Wt(3,1) * (u(4)-u(3))
Csys(2) = (Wt(1,2) * (u(2)-u(1)) + Wt(2,2) * (u(3)-u(2))) + Wt(3,2) * (u(4)-u(3))
Csys(3) = (Wt(1,3) * (u(2)-u(1)) + Wt(2,3) * (u(3)-u(2))) + Wt(3,3) * (u(4)-u(3))
Csys(4) = (Wt(1,4) * (u(2)-u(1)) + Wt(2,4) * (u(3)-u(2))) + Wt(3,4) * (u(4)-u(3))

! endif ! End of non-uniform layer thickness branch.

end subroutine end_value_h4

!> Value of PPM_H4_2019 reconstruction at a point in cell k [A]
real function f(this, k, x)
class(PPM_H4_2019), intent(in) :: this !< This reconstruction
Expand Down Expand Up @@ -310,7 +562,7 @@ logical function unit_tests(this, verbose, stdout, stderr)

end function unit_tests

!> \namespace recon1d_ppm_c4
!> \namespace recon1d_ppm_c4_2019
!!

end module Recon1d_PPM_H4_2019

0 comments on commit 1793985

Please sign in to comment.