Recon1d_EPPM_CWK.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Piecewise Parabolic Method 1D reconstruction in model index space with linear
6!! extrapolation for first and last cells
7!!
8!! This implementation of PPM follows Colella and Woodward, 1984, using uniform thickness
9!! and with cells resorting to PCM for local extrema. First and last cells use a PLM
10!! representation with slope set by matching the edge of the first interior cell.
12
13use recon1d_type, only : recon1d, testing
14use recon1d_ppm_cwk, only : ppm_cwk
15
16implicit none ; private
17
18public eppm_cwk, testing
19
20!> PPM reconstruction in index space (no grid dependence) with linear extrapolation
21!! for first and last cells.
22!!
23!! Implemented by extending recon1d_ppm_cwk.
24!!
25!! The source for the methods ultimately used by this class are:
26!! - init() -> recon1d_ppm_cwk.init()
27!! - reconstruct() *locally defined
28!! - average() -> recon1d_ppm_cwk.average()
29!! - f() -> recon1d_ppm_cwk.f()
30!! - dfdx() -> recon1d_ppm_cwk.dfdx()
31!! - check_reconstruction() -> recon1d_ppm_cwk.check_reconstruction()
32!! - unit_tests() *locally defined
33!! - destroy() -> recon1d_ppm_cwk.destroy()
34!! - remap_to_sub_grid() -> recon1d_type.remap_to_sub_grid()
35!! - init_parent() -> recon1d_ppm_cwk.init()
36!! - reconstruct_parent() -> recon1d_ppm_cwk.reconstruct()
37type, extends (ppm_cwk) :: eppm_cwk
38
39contains
40 !> Implementation of the EPPM_CWK reconstruction
41 procedure :: reconstruct => reconstruct
42 !> Implementation of unit tests for the EPPM_CWK reconstruction
43 procedure :: unit_tests => unit_tests
44
45end type eppm_cwk
46
47contains
48
49!> Calculate a 1D EPPM_CWK reconstructions based on h(:) and u(:)
50subroutine reconstruct(this, h, u)
51 class(eppm_cwk), intent(inout) :: this !< This reconstruction
52 real, intent(in) :: h(*) !< Grid spacing (thickness) [typically H]
53 real, intent(in) :: u(*) !< Cell mean values [A]
54 ! Local variables
55 real :: dul, dur ! Left and right cell PLM slopes [A]
56 real :: u0, u1, u2 ! Far left, left, and right cell values [A]
57 real :: edge ! Edge value between cell k-1 and k [A]
58 real :: u_min, u_max ! Minimum and maximum value across edge [A]
59 real :: a6 ! Colella and Woodward curvature [A]
60 real :: du ! Difference between edges across cell [A]
61 real :: slp(this%n) ! PLM slope [A]
62 real, parameter :: one_sixth = 1. / 6. ! 1/6 [nondim]
63 integer :: k, n
64
65 n = this%n
66
67 call this%reconstruct_parent( h, u )
68
69 ! Extrapolate in first cell
70 this%ur(1) = this%ul(2) ! Assume ur=ul on right edge
71 this%ul(1) = u(1) + ( u(1) - this%ur(1) ) ! Linearly extrapolat across cell
72
73 ! Extrapolate in last cell
74 this%ul(n) = this%ur(n-1) ! Assume ul=ur on left edge
75 this%ur(n) = u(n) + ( u(n) - this%ul(n) ) ! Linearly extrapolat across cell
76
77end subroutine reconstruct
78
79!> Runs EPPM_CWK reconstruction unit tests and returns True for any fails, False otherwise
80logical function unit_tests(this, verbose, stdout, stderr)
81 class(eppm_cwk), intent(inout) :: this !< This reconstruction
82 logical, intent(in) :: verbose !< True, if verbose
83 integer, intent(in) :: stdout !< I/O channel for stdout
84 integer, intent(in) :: stderr !< I/O channel for stderr
85 ! Local variables
86 real, allocatable :: ul(:), ur(:), um(:) ! test values [A]
87 real, allocatable :: ull(:), urr(:) ! test values [A]
88 type(testing) :: test ! convenience functions
89 integer :: k
90
91 call test%set( stdout=stdout ) ! Sets the stdout channel in test
92 call test%set( stderr=stderr ) ! Sets the stderr channel in test
93 call test%set( verbose=verbose ) ! Sets the verbosity flag in test
94
95 if (verbose) write(stdout,'(a)') 'EPPM_CWK:unit_tests testing with linear fn'
96
97 call this%init(5)
98 call test%test( this%n /= 5, 'Setting number of levels')
99 allocate( um(5), ul(5), ur(5), ull(5), urr(5) )
100
101 ! Straight line, f(x) = x , or f(K) = 2*K
102 call this%reconstruct( (/2.,2.,2.,2.,2./), (/1.,4.,7.,10.,13./) )
103 call test%real_arr(5, this%u_mean, (/1.,4.,7.,10.,13./), 'Setting cell values')
104 call test%real_arr(5, this%ul, (/-0.5,2.5,5.5,8.5,11.5/), 'Left edge values')
105 call test%real_arr(5, this%ur, (/2.5,5.5,8.5,11.5,14.5/), 'Right edge values')
106
107 do k = 1, 5
108 ul(k) = this%f(k, 0.)
109 um(k) = this%f(k, 0.5)
110 ur(k) = this%f(k, 1.)
111 enddo
112 call test%real_arr(5, ul, this%ul, 'Evaluation on left edge')
113 call test%real_arr(5, um, (/1.,4.,7.,10.,13./), 'Evaluation in center')
114 call test%real_arr(5, ur, this%ur, 'Evaluation on right edge')
115
116 do k = 1, 5
117 ul(k) = this%dfdx(k, 0.)
118 um(k) = this%dfdx(k, 0.5)
119 ur(k) = this%dfdx(k, 1.)
120 enddo
121 ! Most of these values are affected by the PLM boundary cells
122 call test%real_arr(5, ul, (/3.,3.,3.,3.,3./), 'dfdx on left edge')
123 call test%real_arr(5, um, (/3.,3.,3.,3.,3./), 'dfdx in center')
124 call test%real_arr(5, ur, (/3.,3.,3.,3.,3./), 'dfdx on right edge')
125
126 do k = 1, 5
127 um(k) = this%average(k, 0.5, 0.75) ! Average from x=0.25 to 0.75 in each cell
128 enddo
129 ! Most of these values are affected by the PLM boundary cells
130 call test%real_arr(5, um, (/1.375,4.375,7.375,10.375,13.375/), 'Return interval average')
131
132 if (verbose) write(stdout,'(a)') 'EPPM_CWK:unit_tests testing with parabola'
133
134 ! x = 2 i i=0 at origin
135 ! f(x) = 3/4 x^2 = (2 i)^2
136 ! f[i] = 3/4 ( 2 i - 1 )^2 on centers
137 ! f[I] = 3/4 ( 2 I )^2 on edges
138 ! f[i] = 1/8 [ x^3 ] for means
139 ! edges: 0, 1, 12, 27, 48, 75
140 ! means: 1, 7, 19, 37, 61
141 ! centers: 0.75, 6.75, 18.75, 36.75, 60.75
142 call this%reconstruct( (/2.,2.,2.,2.,2./), (/1.,7.,19.,37.,61./) )
143 do k = 1, 5
144 ul(k) = this%f(k, 0.)
145 um(k) = this%f(k, 0.5)
146 ur(k) = this%f(k, 1.)
147 enddo
148 call test%real_arr(5, ul, (/-1.,3.,12.,27.,48./), 'Return left edge')
149 call test%real_arr(5, um, (/1.,6.75,18.75,36.75,61./), 'Return center')
150 call test%real_arr(5, ur, (/3.,12.,27.,48.,74./), 'Return right edge')
151
152 ! x = 3 i i=0 at origin
153 ! f(x) = x^2 / 3 = 3 i^2
154 ! f[i] = [ ( 3 i )^3 - ( 3 i - 3 )^3 ] i=1,2,3,4,5
155 ! means: 1, 7, 19, 37, 61
156 ! edges: 0, 3, 12, 27, 48, 75
157 call this%reconstruct( (/3.,3.,3.,3.,3./), (/1.,7.,19.,37.,61./) )
158 do k = 1, 5
159 ul(k) = this%f(k, 0.)
160 um(k) = this%f(k, 0.5)
161 ur(k) = this%f(k, 1.)
162 enddo
163 call test%real_arr(5, ul, (/-1.,3.,12.,27.,48./), 'Return left edge')
164 call test%real_arr(5, um, (/1.,6.75,18.75,36.75,61./), 'Return center')
165 call test%real_arr(5, ur, (/3.,12.,27.,48.,74./), 'Return right edge')
166
167 call this%destroy()
168 deallocate( um, ul, ur, ull, urr )
169
170 unit_tests = test%summarize('EPPM_CWK:unit_tests')
171
172end function unit_tests
173
174!> \namespace recon1d_eppm_cwk
175!!
176
177end module recon1d_eppm_cwk