diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index 8b3bbacdce..297e45b227 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -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_CW, only : PPM_CW use Recon1d_PPM_H4_2019, only : PPM_H4_2019 implicit none ; private @@ -1760,6 +1761,9 @@ subroutine setReconstructionType(string,CS) case ("C_EMPLM_WA_POLY") allocate( EMPLM_WA_poly :: CS%reconstruction ) CS%remapping_scheme = REMAPPING_VIA_CLASS + case ("C_PPM_CW") + allocate( PPM_CW :: CS%reconstruction ) + CS%remapping_scheme = REMAPPING_VIA_CLASS case ("C_PPM_H4_2019") allocate( PPM_H4_2019 :: CS%reconstruction ) CS%remapping_scheme = REMAPPING_VIA_CLASS @@ -1820,6 +1824,7 @@ logical function remapping_unit_tests(verbose, num_comp_samp) type(MPLM_WA_poly) :: MPLM_WA_poly type(EMPLM_WA_poly) :: EMPLM_WA_poly type(PPM_H4_2019) :: PPM_H4_2019 + type(PPM_CW) :: PPM_CW call test%set( verbose=verbose ) ! Sets the verbosity flag in test call test%set( stop_instantly=.true. ) ! While debugging @@ -2438,7 +2443,9 @@ 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_CW%unit_tests(verbose, test%stdout, test%stderr), 'PPM_CW unit test') call test%test( PPM_H4_2019%unit_tests(verbose, test%stdout, test%stderr), 'PPM_H4_2019 unit test') +stop ! Randomized, brute force tests ntests = 3000 diff --git a/src/ALE/Recon1d_PPM_CW.F90 b/src/ALE/Recon1d_PPM_CW.F90 new file mode 100644 index 0000000000..b2d5998d16 --- /dev/null +++ b/src/ALE/Recon1d_PPM_CW.F90 @@ -0,0 +1,407 @@ +!> Piecewise Parabolic Method 1D reconstruction following Colella and Woodward, 1984 +!! +!! This implementation of PPM follows Colella and Woodward, 1984, with cells resorting to PCM for +!! extrema including first and last cells in column. +module Recon1d_PPM_CW + +! This file is part of MOM6. See LICENSE.md for the license. + +use Recon1d_type, only : Recon1d, testing +use Recon1d_PLM_CW, only : PLM_CW + +implicit none ; private + +public PPM_CW, 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 +!! average() *locally defined +!! f() *locally defined +!! dfdx() *locally defined +!! check_reconstruction() *locally defined +!! unit_tests() *locally defined +!! destroy() *locally defined +!! remap_to_sub_grid() -> Recon1d%remap_to_sub_grid() +!! init_parent() -> init() +!! reconstruct_parent() -> reconstruct() +type, extends (Recon1d) :: PPM_CW + + real, allocatable :: ul(:) !< Left edge value [A] + real, allocatable :: ur(:) !< Right edge value [A] + type(PLM_CW) :: PLM !< The PLM reconstruction used to estimate edge values + +contains + !> Implementation of the PPM_CW initialization + procedure :: init => init + !> Implementation of the PPM_CW reconstruction + procedure :: reconstruct => reconstruct + !> Implementation of the PPM_CW average over an interval [A] + procedure :: average => average + !> Implementation of evaluating the PPM_CW reconstruction at a point [A] + procedure :: f => f + !> Implementation of the derivative of the PPM_CW reconstruction at a point [A] + procedure :: dfdx => dfdx + !> Implementation of deallocation for PPM_CW + procedure :: destroy => destroy + !> Implementation of check reconstruction for the PPM_CW reconstruction + procedure :: check_reconstruction => check_reconstruction + !> Implementation of unit tests for the PPM_CW 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_CW + +contains + +!> Initialize a 1D PPM_CW reconstruction for n cells +subroutine init(this, n, h_neglect, check) + class(PPM_CW), 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 incurs an extra store of u_mean but by using PCM_CW + ! we avoid duplicating and testing more code + call this%PLM%init( n, h_neglect=h_neglect, check=check ) + + 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_CW reconstructions based on h(:) and u(:) +subroutine reconstruct(this, h, u) + class(PPM_CW), 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 :: h0, h1, h2, h3 ! Cell thickness h(k-2), h(k-1), h(k), h(k+1) in K loop [H] + real :: d12 ! h1 + h2 but used in the denominator so include h_neglect [H] + real :: h01_h112, h23_h122 ! Approximately 2/3 [nondim] + real :: ddh ! Approximately 0 [nondim] + real :: I_h12, I_h0123 ! Reciprocals of d12 and sum(h) [H-1] + real :: dul, dur ! Left and right cell PLM slopes [A] + real :: u0, u1, u2 ! Far left, left, and right cell values [A] + real :: edge ! Edge value between cell k-1 and k [A] + real :: u_min, u_max ! Minimum and maximum value across edge [A] + real :: a6 ! Colella and Woodward curvature [A] + real :: du ! Difference between edges across cell [A] + real :: slp(this%n) ! PLM slope [A] + integer :: k, n + + n = this%n + +print *,'u=',u(1:n) + ! First populate the PLM reconstructions + call this%PLM%reconstruct( h, u ) + do k = 1, n + slp(k) = this%PLM%ur(k) - this%PLM%ul(k) + enddo + slp(1) = 2.0 * ( this%PLM%ul(2) - u(1) ) + slp(n) = 2.0 * ( u(n) - this%PLM%ur(n-1) ) + + do K = 2, n ! K=2 is interface between cells 1 and 2 + h0 = h( max( 1, k-2 ) ) ! This treatment implies a virtual mirror cell at k=0 + h1 = h(k-1) + h2 = h(k) + h3 = h( min( n, k+1 ) ) ! This treatment implies a virtual mirror cell at k=n+1 + d12 = ( h1 + h2 ) + this%h_neglect ! d12 is only ever used in the denominator + h01_h112 = ( h0 + h1 ) / ( h1 + d12 ) ! When uniform -> 2/3 + h23_h122 = ( h2 + h3 ) / ( d12 + h2 ) ! When uniform -> 2/3 + ddh = h01_h112 - h23_h122 ! When uniform -> 0 + I_h12 = 1.0 / d12 ! When uniform -> 1/(2h) + I_h0123 = 1.0 / ( d12 + ( h0 + h3 ) ) ! When uniform -> 1/(4h) + dul = slp(k-1) + dur = slp(k) + u2 = u(k) + u1 = u(k-1) + print *,k,'uk-1=',u1,'du-=',dul,'du+=',dur,'uk=',u2 + edge = I_h12 * ( h2 * u1 + h1 * u2 ) & ! 1/2 u1 + 1/2 u2 + + I_h0123 * ( 2.0 * h1 * h2 * I_h12 * ddh * ( u2 - u1 ) & ! 0 + + ( h2 * h23_h122 * dul - h1 * h01_h112 * dur ) ) ! 1/12 dul - 1/12 dur + u_min = min( u1, u2 ) + u_max = max( u1, u2 ) + edge = max( min( edge, u_max), u_min ) ! Unclear if we need this bounding in the interior + this%ur(k-1) = edge + this%ul(k) = edge +print *,k,'edge=',edge,'mid=',I_h12 * ( h2 * u1 + h1 * u2 ) +print *,k,' O(4) corr=', I_h0123 * ( h2 * h23_h122 * dul - h1 * h01_h112 * dur ) +!if (k>2 .and. k du * du ) then ! Extrema on right + edge = 3.0 * u1 - 2.0 * this%ur(k) + ! u_min = min( u0, u1 ) + ! u_max = max( u0, u1 ) + ! edge = max( min( edge, u_max), u_min ) + this%ul(k) = edge + elseif ( du * a6 < - du * du ) then ! Extrema on left + edge = 3.0 * u1 - 2.0 * this%ul(k) + ! u_min = min( u1, u2 ) + ! u_max = max( u1, u2 ) + ! edge = max( min( edge, u_max), u_min ) + this%ur(k) = edge + endif + enddo + + ! After the limiter, are ur and ul bounded???? -AJA + + ! Store mean + do k = 1, n + this%u_mean(k) = u(k) + enddo + +end subroutine reconstruct + +!> Value of PPM_CW reconstruction at a point in cell k [A] +real function f(this, k, x) + class(PPM_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 :: a6 ! Collela and Woordward curvature parameter [A] + real :: u_a, u_b ! Two estimate of f [A] + real :: lmx ! 1 - x [nondim] + real :: wb ! Weight based on x [nondim] + + 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 ) ) + lmx = 1.0 - xc + + ! This expression for u_a can overshoot u_r but is good for x<<1 + 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, xc - 0.5 ) ! = 1 @ x=0, = 0 @ x=1 + f = ( ( 1. - wb ) * u_a ) + ( wb * u_b ) + +end function f + +!> Derivative of PPM_CW reconstruction at a point in cell k [A] +real function dfdx(this, k, x) + class(PPM_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 :: 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 * xc - 1.0 ) + +end function dfdx + +!> Average between xa and xb for cell k of a 1D PPM reconstruction [A] +real function average(this, k, xa, xb) + class(PPM_CW), 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_CW reconstruction +subroutine destroy(this) + class(PPM_CW), intent(inout) :: this !< This reconstruction + + deallocate( this%u_mean, this%ul, this%ur ) + +end subroutine destroy + +!> Checks the PPM_CW reconstruction for consistency +logical function check_reconstruction(this, h, u) + class(PPM_CW), 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_CW reconstruction unit tests and returns True for any fails, False otherwise +logical function unit_tests(this, verbose, stdout, stderr) + class(PPM_CW), 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 +call test%set( stop_instantly=.true. ) ! debugging + + if (verbose) write(stdout,'(a)') 'PPM_CW: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) ) + + ! Straight line, f(x) = x , or f(K) = 2*K + call this%reconstruct( (/2.,2.,2.,2.,2./), (/1.,4.,7.,10.,13./) ) + call test%real_arr(5, this%u_mean, (/1.,4.,7.,10.,13./), 'Setting cell values') + ! Without PLM extrapolation we get l(2)=2 and r(4)=12 due to PLM=0 in boundary cells. -AJA + call test%real_arr(5, this%ul, (/1.,2.5,5.5,8.5,13./), 'Left edge values') + call test%real_arr(5, this%ur, (/1.,5.5,8.5,11.5,13./), 'Right edge values') + + do k = 1, 5 + ul(k) = this%f(k, 0.) + um(k) = this%f(k, 0.5) + ur(k) = this%f(k, 1.) + enddo + call test%real_arr(5, ul, this%ul, 'Evaluation on left edge') + call test%real_arr(5, um, (/1.,4.,7.,10.,13./), 'Evaluation in center') + call test%real_arr(5, ur, this%ur, 'Evaluation on right edge') + + do k = 1, 5 + ul(k) = this%dfdx(k, 0.) + um(k) = this%dfdx(k, 0.5) + ur(k) = this%dfdx(k, 1.) + enddo + ! Most of these values are affected by the PLM boundary cells + call test%real_arr(5, ul, (/0.,3.,3.,3.,0./), 'dfdx on left edge') + call test%real_arr(5, um, (/0.,3.,3.,3.,0./), 'dfdx in center') + call test%real_arr(5, ur, (/0.,3.,3.,3.,0./), 'dfdx on 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 + ! Most of these values are affected by the PLM boundary cells + call test%real_arr(5, um, (/1.,4.375,7.375,10.375,13./), 'Return interval average') + + if (verbose) write(stdout,'(a)') 'PPM_CW:unit_tests testing with parabola' + + ! x = 2 i i=0 at origin + ! f(x) = x^2 = (2 i)^2 + ! f[i] = ( 2 i - 1 )^2 on centers + ! f[I] = ( 2 I )^2 on edges + ! means: 1, 9, 25, 49, 81 + ! edges: 0, 4, 16, 36, 64, 100 + ! slopes 0, 12, 20, 28, 0 + ! mid-points: 5, 17, 37, 65 + ! (dur-dul)/12: -1,-2/3,-2/3,28/12 + ! edge: 3, +print *,'#########################################################################################' + call this%reconstruct( (/2.,2.,2.,2.,2./), (/1.,9.,25.,49.,81./) ) + do k = 1, 5 + ul(k) = this%f(k, 0.) + um(k) = this%f(k, 0.5) + ur(k) = this%f(k, 1.) + 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) + + ! 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 +! ul(k) = this%f(k, 0.) +! um(k) = this%f(k, 0.5) +! ur(k) = this%f(k, 1.) +! 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) + + call this%destroy() + deallocate( um, ul, ur, ull, urr ) + + unit_tests = test%summarize('PPM_CW:unit_tests') + +end function unit_tests + +!> \namespace recon1d_ppm_c4 +!! + +end module Recon1d_PPM_CW