Skip to content

Commit

Permalink
Added draft of PPM_H4_2019
Browse files Browse the repository at this point in the history
  • Loading branch information
@adcroft
adcroft committed Aug 8, 2024
1 parent 3f8d067 commit 9984bf8
Show file tree
Hide file tree
Showing 2 changed files with 349 additions and 0 deletions.
3 changes: 3 additions & 0 deletions src/ALE/MOM_remapping.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@ module MOM_remapping
use Recon1d_EMPLM_WA, only : EMPLM_WA
use Recon1d_MPLM_WA_poly, only : MPLM_WA_poly
use Recon1d_EMPLM_WA_poly, only : EMPLM_WA_poly
use Recon1d_PPM_H4_2019, only : PPM_H4_2019

implicit none ; private

Expand Down Expand Up @@ -1719,6 +1720,7 @@ logical function remapping_unit_tests(verbose, num_comp_samp)
type(EMPLM_WA) :: EMPLM_WA
type(MPLM_WA_poly) :: MPLM_WA_poly
type(EMPLM_WA_poly) :: EMPLM_WA_poly
type(PPM_H4_2019) :: PPM_H4_2019

call test%set( verbose=verbose ) ! Sets the verbosity flag in test

Expand Down Expand Up @@ -2336,6 +2338,7 @@ logical function remapping_unit_tests(verbose, num_comp_samp)
call test%test( MPLM_WA_poly%unit_tests(verbose, test%stdout, test%stderr), 'MPLM_WA_poly unit test')
call test%test( EMPLM_WA_poly%unit_tests(verbose, test%stdout, test%stderr), 'EMPLM_WA_poly unit test')
call test%test( PLM_CW%unit_tests(verbose, test%stdout, test%stderr), 'PLM_CW unit test')
call test%test( PPM_H4_2019%unit_tests(verbose, test%stdout, test%stderr), 'PPM_H4_2019 unit test')

n0 = 8
allocate( h0(n0), u0(n0) )
Expand Down
346 changes: 346 additions & 0 deletions src/ALE/Recon1d_PPM_H4_2019.F90
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,346 @@
!> Piecewise Parabolic Method 1D reconstruction with h4 interpolation for edges
!!
!! 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.
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

public PPM_H4_2019, testing

!> PPM reconstruction following White and Adcroft, 2008
!!
!! The source for the methods ultimately used by this class are:
!! init() *locally defined
!! reconstruct() *locally defined
!! lr_edge() *locally defined
!! average() *locally defined
!! inv_f() *locally defined
!! check_reconstruction() *locally defined
!! unit_tests() *locally defined
!! destroy() *locally defined
!! cell_mean() -> Recon1d%cell_mean()
!! remap_to_sub_grid() -> Recon1d%remap_to_sub_grid()
!! init_parent() -> init()
!! reconstruct_parent() -> reconstruct()
type, extends (Recon1d) :: PPM_H4_2019

real, allocatable :: ul(:) !< Left edge value [A]
real, allocatable :: ur(:) !< Right edge value [A]

contains
!> Implementation of the PPM_H4_2019 initialization
procedure :: init => init
!> Implementation of the PPM_H4_2019 reconstruction
procedure :: reconstruct => reconstruct
!> Implementation of function returning the PPM_H4_2019 edge values
procedure :: lr_edge => lr_edge
!> Implementation of the PPM_H4_2019 average over an interval [A]
procedure :: average => average
!> Implementation of finding the PPM_H4_2019 position of a value
procedure :: inv_f => inv_f
!> Implementation of deallocation for PPM_H4_2019
procedure :: destroy => destroy
!> Implementation of check reconstruction for the PPM_H4_2019 reconstruction
procedure :: check_reconstruction => check_reconstruction
!> Implementation of unit tests for the PPM_H4_2019 reconstruction
procedure :: unit_tests => unit_tests

!> Duplicate interface to init()
procedure :: init_parent => init
!> Duplicate interface to reconstruct()
procedure :: reconstruct_parent => reconstruct

end type PPM_H4_2019

contains

!> Initialize a 1D PPM_H4_2019 reconstruction for n cells
subroutine init(this, n, h_neglect, check)
class(PPM_H4_2019), intent(out) :: this !< This reconstruction
integer, intent(in) :: n !< Number of cells in this column
real, optional, intent(in) :: h_neglect !< A negligibly small width used in cell reconstructions [H]
logical, optional, intent(in) :: check !< If true, enable some consistency checking

this%n = n

allocate( this%u_mean(n) )
allocate( this%ul(n) )
allocate( this%ur(n) )

this%h_neglect = tiny( this%u_mean(1) )
if (present(h_neglect)) this%h_neglect = h_neglect
this%check = .false.
if (present(check)) this%check = check

end subroutine init

!> Calculate a 1D PPM_H4_2019 reconstructions based on h(:) and u(:)
subroutine reconstruct(this, h, u)
class(PPM_H4_2019), intent(inout) :: this !< This reconstruction
real, intent(in) :: h(*) !< Grid spacing (thickness) [typically H]
real, intent(in) :: u(*) !< Cell mean values [A]
! Local variables
real :: slp ! The PLM slopes (difference across cell) [A]
real :: sigma_l, sigma_c, sigma_r ! Left, central and right slope estimates as
! differences across the cell [A]
real :: u_min, u_max ! Minimum and maximum value across cell [A]
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 :: edge_values(this%n,2) ! Edge values [A]
real :: ppoly_coef(this%n,3) ! Polynomial coefficients [A]

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 )

! Store reconstruction
do k = 1, n
this%u_mean(k) = u(k)
this%ul(k) = edge_values(k,1)
this%ur(k) = edge_values(k,2)
enddo

end subroutine reconstruct

!> Returns the left/right edge values in cell k of a 1D PPM reconstruction
subroutine lr_edge(this, k, u_left, u_right)
class(PPM_H4_2019), intent(in) :: this !< This reconstruction
integer, intent(in) :: k !< Cell number
real, intent(out) :: u_left !< Left edge value [A]
real, intent(out) :: u_right !< Right edge value [A]

u_left = this%ul(k)
u_right = this%ur(k)

end subroutine lr_edge

!> Position at which 1D PPM_H4_2019 reconstruction has a value f in cell k [0..1]
real function inv_f(this, k, f)
class(PPM_H4_2019), intent(in) :: this !< This reconstruction
integer, intent(in) :: k !< Cell number
real, intent(in) :: f !< Value of reconstruction to solve for [A]
real :: dfdx_sign ! sign of dfdx, i.e. -1 or +1 [nondim]
real :: sf ! Values with potentially flipped sign [A]
real :: du ! Difference across cell [A]
real :: a6 ! The Colella and Woodward PPM curviture parameter [A]
real :: qA, qB, qC ! Respective multiplier of x^2, x, 1 in quadratic [nondim]
real :: discr ! Discriminant in quadratic problem [nondim]

du = this%ur(k) - this%ul(k)

! This is the sign of dfdx across cell k
if (abs(du) <= 0.) then
dfdx_sign = sign( 1., this%u_mean(min(this%n,k+1)) - this%u_mean(max(1,k-1)) )
else
dfdx_sign = sign( 1., du )
endif
sf = f * dfdx_sign

inv_f = 0.5 ! Avoid compiler warnings

if ( sf < this%ul(k) * dfdx_sign ) then
inv_f = 0.
elseif ( sf > this%ur(k) * dfdx_sign ) then
inv_f = 1.
elseif ( du == 0. ) then
inv_f = 0.5 ! Same as for PCM
else ! ul < f < ur or ul > f > ur
! Invert PPM since we know f is between ul and ur
a6 = 3.0 * ( ( this%u_mean(k) - this%ul(k) ) + ( this%u_mean(k) - this%ur(k) ) )
if ( abs(a6) > 0. ) then ! u_mean not equal to (ul+ur)/2
qA = - a6
qB = du + a6
qC = this%ul(k) - f
discr = sqrt( (qB * qB) - 4. * ( qA * qC ) )
if ( qB >= 0. ) then
inv_f = - 0.5 * ( qB + discr ) / qA
inv_f = - 2.0 * qC / ( qB + discr )
else ! qB < 0
inv_f = - 0.5 * ( qB + discr ) / qA
inv_f = 2.0 * qC / ( - qB + discr )
endif
else ! a6==0
inv_f = ( f - this%ul(k) ) / du
endif
endif

end function inv_f

!> Average between xa and xb for cell k of a 1D PPM reconstruction [A]
real function average(this, k, xa, xb)
class(PPM_H4_2019), intent(in) :: this !< This reconstruction
integer, intent(in) :: k !< Cell number
real, intent(in) :: xa !< Start of averaging interval on element (0 to 1)
real, intent(in) :: xb !< End of averaging interval on element (0 to 1)
real :: xapxb ! A sum of fracional positions [nondim]
real :: mx, Ya, Yb, my ! Various fractional positions [nondim]
real :: u_a, u_b ! Values at xa and xb [A]
real :: xa2pxb2, xa2b2ab, Ya2b2ab ! Sums of squared fractional positions [nondim]
real :: a_L, a_R, u_c, a_c ! Values of the polynomial at various locations [A]

mx = 0.5 * ( xa + xb )
a_L = this%ul(k)
a_R = this%ur(k)
u_c = this%u_mean(k)
a_c = 0.5 * ( ( u_c - a_L ) + ( u_c - a_R ) ) ! a_6 / 6
if (mx<0.5) then
! This integration of the PPM reconstruction is expressed in distances from the left edge
xa2b2ab = (xa*xa+xb*xb)+xa*xb
average = a_L + ( ( a_R - a_L ) * mx &
+ a_c * ( 3. * ( xb + xa ) - 2.*xa2b2ab ) )
else
! This integration of the PPM reconstruction is expressed in distances from the right edge
Ya = 1. - xa
Yb = 1. - xb
my = 0.5 * ( Ya + Yb )
Ya2b2ab = (Ya*Ya+Yb*Yb)+Ya*Yb
average = a_R + ( ( a_L - a_R ) * my &
+ a_c * ( 3. * ( Yb + Ya ) - 2.*Ya2b2ab ) )
endif

end function average

!> Deallocate the PPM_H4_2019 reconstruction
subroutine destroy(this)
class(PPM_H4_2019), intent(inout) :: this !< This reconstruction

deallocate( this%u_mean, this%ul, this%ur )

end subroutine destroy

!> Checks the PPM_H4_2019 reconstruction for consistency
logical function check_reconstruction(this, h, u)
class(PPM_H4_2019), intent(in) :: this !< This reconstruction
real, intent(in) :: h(*) !< Grid spacing (thickness) [typically H]
real, intent(in) :: u(*) !< Cell mean values [A]
! Local variables
integer :: k

check_reconstruction = .false.

do k = 1, this%n
if ( abs( this%u_mean(k) - u(k) ) > 0. ) check_reconstruction = .true.
enddo

! Check implied curvature
do k = 1, this%n
if ( ( this%u_mean(k) - this%ul(k) ) * ( this%ur(k) - this%u_mean(k) ) < 0. ) check_reconstruction = .true.
enddo

! Check bounding of right edges
do K = 1, this%n-1
if ( ( this%ur(k) - this%u_mean(k) ) * ( this%u_mean(k+1) - this%ur(k) ) < 0. ) check_reconstruction = .true.
enddo

! Check bounding of left edges
do K = 2, this%n
if ( ( this%u_mean(k) - this%ul(k) ) * ( this%ul(k) - this%u_mean(k-1) ) < 0. ) check_reconstruction = .true.
enddo

end function check_reconstruction

!> Runs PPM_H4_2019 reconstruction unit tests and returns True for any fails, False otherwise
logical function unit_tests(this, verbose, stdout, stderr)
class(PPM_H4_2019), intent(inout) :: this !< This reconstruction
logical, intent(in) :: verbose !< True, if verbose
integer, intent(in) :: stdout !< I/O channel for stdout
integer, intent(in) :: stderr !< I/O channel for stderr
! Local variables
real, allocatable :: ul(:), ur(:), um(:) ! test values [A]
real, allocatable :: ull(:), urr(:) ! test values [A]
type(testing) :: test ! convenience functions
integer :: k

call test%set( stdout=stdout ) ! Sets the stdout channel in test
call test%set( stderr=stderr ) ! Sets the stderr channel in test
call test%set( verbose=verbose ) ! Sets the verbosity flag in test

if (verbose) write(stdout,'(a)') 'PPM_H4_2019:unit_tests testing with linear fn'

call this%init(5)
call test%test( this%n /= 5, 'Setting number of levels')
allocate( um(5), ul(5), ur(5), ull(5), urr(5) )

call this%reconstruct( (/2.,2.,2.,2.,2./), (/1.,3.,5.,7.,9./) )
call test%real_arr(5, this%u_mean, (/1.,3.,5.,7.,9./), 'Setting cell values')
do k = 1, 5
um(k) = this%cell_mean(k)
enddo
call test%real_arr(5, um, (/1.,3.,5.,7.,9./), 'Return cell mean')

do k = 1, 5
call this%lr_edge(k, ul(k), ur(k))
enddo
! Values of 2 at k=2 differs by roundoff
call test%real_arr(5, ul, (/1.,2.,4.,6.,9./), 'Return left edge', robits=2)
call test%real_arr(5, ur, (/1.,4.,6.,8.,9./), 'Return right edge')

do k = 1, 5
um(k) = this%average(k, 0.5, 0.75) ! Average from x=0.25 to 0.75 in each cell
enddo
call test%real_arr(5, um, (/1.,3.25,5.25,7.25,9./), 'Return interval average')

do k = 1, 5
ull(k) = this%inv_f(k, real(2*k-1)-1.0)
ul(k) = this%inv_f(k, real(2*k-1)-0.5)
um(k) = this%inv_f(k, real(2*k-1))
ur(k) = this%inv_f(k, real(2*k-1)+0.5)
urr(k) = this%inv_f(k, real(2*k-1)+1.0)
enddo
call test%real_arr(5, ull, (/0.,0.,0.,0.,0./), 'Return position of f<<', tol=4e-16)
call test%real_arr(5, ul, (/0.,0.25,0.25,0.25,0./), 'Return position of f<', robits=4)
call test%real_arr(5, um, (/0.5,0.5,0.5,0.5,0.5/), 'Return position of f=', robits=1)
call test%real_arr(5, ur, (/1.,0.75,0.75,0.75,1./), 'Return position of f>', robits=1)
call test%real_arr(5, urr, (/1.,1.,1.,1.,1./), 'Return position of f>>')

if (verbose) write(stdout,'(a)') 'PPM_H4_2019:unit_tests testing with parabola'

! x = 3 i i=0 at origin
! f(x) = x^2 / 3 = 3 i^2
! f[i] = [ ( 3 i )^3 - ( 3 i - 3 )^3 ] i=1,2,3,4,5
! means: 1, 7, 19, 37, 61
! edges: 0, 3, 12, 27, 48, 75
call this%reconstruct( (/3.,3.,3.,3.,3./), (/1.,7.,19.,37.,61./) )
do k = 1, 5
call this%lr_edge(k, ul(k), ur(k))
enddo
call test%real_arr(5, ul, (/1.,3.,12.,27.,61./), 'Return left edge', robits=2)
call test%real_arr(5, ur, (/1.,12.,27.,48.,61./), 'Return right edge', robits=1)

! do k = 1, 5
! ull(k) = this%inv_f(k, 3.*(real(k)-1.0)**2 )
! ul(k) = this%inv_f(k, 3.*(real(k)-0.75)**2 )
! um(k) = this%inv_f(k, 3.*(real(k)-0.5)**2 )
! ur(k) = this%inv_f(k, (real(k)-0.5)**2 )
! urr(k) = this%inv_f(k, (real(k)-0.5)**2 )
! enddo
! call test%real_arr(5, ull, (/0.,0.,0.,0.,0./), 'Return position of f<<', tol=5e-16)
! call test%real_arr(5, ul, (/0.,0.25,0.25,0.25,0./), 'Return position of f<', robits=2)
! call test%real_arr(5, um, (/0.5,0.5,0.5,0.5,0.5/), 'Return position of f=', robits=1)
! call test%real_arr(5, ur, (/1.,0.75,0.75,0.75,1./), 'Return position of f>', robits=1)
! call test%real_arr(5, urr, (/1.,1.,1.,1.,1./), 'Return position of f>>')

call this%destroy()
deallocate( um, ul, ur, ull, urr )

unit_tests = test%summarize('PPM_H4_2019:unit_tests')

end function unit_tests

!> \namespace recon1d_ppm_c4
!!

end module Recon1d_PPM_H4_2019

0 comments on commit 9984bf8

Please sign in to comment.