Recon1d_PPM_H4_2018.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 with h4 interpolation for edges (2018 version)
6!!
7!! This implementation of PPM follows White and Adcroft 2008 \cite white2008, with cells
8!! resorting to PCM for extrema including first and last cells in column.
9!! This scheme differs from Colella and Woodward, 1984 \cite colella1984, in the method
10!! of first estimating the fourth-order accurate edge values.
11!! This uses numerical expressions that predate a 2019 refactoring.
12!! The first and last cells are always limited to PCM.
14
15use recon1d_ppm_h4_2019, only : ppm_h4_2019, testing
18
19implicit none ; private
20
21public ppm_h4_2018, testing
22
23!> PPM reconstruction following White and Adcroft, 2008
24!!
25!! Implemented by extending recon1d_ppm_h4_2019.
26!!
27!! The source for the methods ultimately used by this class are:
28!! - init() -> recon1d_ppm_h4_2019.init()
29!! - reconstruct() *locally defined
30!! - average() -> recon1d_ppm_h4_2019.average()
31!! - f() -> recon1d_ppm_h4_2019.f()
32!! - dfdx() -> recon1d_ppm_h4_2019.dfdx()
33!! - check_reconstruction() -> recon1d_ppm_h4_2019.check_reconstruction()
34!! - unit_tests() *locally defined
35!! - destroy() -> recon1d_ppm_h4_2019.destroy()
36!! - remap_to_sub_grid() -> recon1d_type.remap_to_sub_grid()
37!! - init_parent() -> recon1d_ppm_h4_2019.init()
38!! - reconstruct_parent() -> recon1d_ppm_h4_2019.reconstruct()
39type, extends (ppm_h4_2019) :: ppm_h4_2018
40
41contains
42 !> Implementation of the PPM_H4_2018 reconstruction
43 procedure :: reconstruct => reconstruct
44 !> Implementation of unit tests for the PPM_H4_2018 reconstruction
45 procedure :: unit_tests => unit_tests
46
47end type ppm_h4_2018
48
49contains
50
51!> Calculate a 1D PPM_H4_2018 reconstructions based on h(:) and u(:)
52subroutine reconstruct(this, h, u)
53 class(ppm_h4_2018), intent(inout) :: this !< This reconstruction
54 real, intent(in) :: h(*) !< Grid spacing (thickness) [typically H]
55 real, intent(in) :: u(*) !< Cell mean values [A]
56 ! Local variables
57 real :: slp ! The PLM slopes (difference across cell) [A]
58 real :: u_l, u_r, u_c ! Left, right, and center values [A]
59 real :: h_l, h_c, h_r ! Thickness of left, center and right cells [H]
60 real :: h0, h1, h2, h3 ! temporary thicknesses [H]
61 real :: h_min ! A minimal cell width [H]
62 real :: f1 ! An auxiliary variable [H]
63 real :: f2 ! An auxiliary variable [A H]
64 real :: f3 ! An auxiliary variable [H-1]
65 real :: et1, et2, et3 ! terms the expression for edge values [A H]
66 real :: dx ! Difference of successive values of x [H]
67 real :: f ! value of polynomial at x in arbitrary units [A]
68 real :: edge_l, edge_r ! Edge values (left and right) [A]
69 real :: expr1, expr2 ! Temporary expressions [A2]
70 real, parameter :: hMinFrac = 1.e-5 !< A minimum fraction for min(h)/sum(h) [nondim]
71 real, dimension(5) :: x ! Coordinate system with 0 at edges [H]
72 real :: edge_values(this%n,2) ! Edge values [A]
73 real :: ppoly_coef(this%n,3) ! Polynomial coefficients [A]
74 real, dimension(4,4) :: A ! Differences in successive positions raised to various powers,
75 ! in units that vary with the second (j) index as [H^j]
76 real, dimension(4) :: B ! The right hand side of the system to solve for C [A H]
77 real, dimension(4) :: C ! The coefficients of a fit polynomial in units that vary
78 ! with the index (j) as [A H^(j-1)]
79 integer :: k, n, j
80
81 n = this%n
82
83 ! Loop on interior cells
84 do k = 3, n-1
85
86 h0 = h(k-2)
87 h1 = h(k-1)
88 h2 = h(k)
89 h3 = h(k+1)
90
91 ! Avoid singularities when consecutive pairs of h vanish
92 if (h0+h1==0.0 .or. h1+h2==0.0 .or. h2+h3==0.0) then
93 h_min = hminfrac*max( this%h_neglect, h0+h1+h2+h3 )
94 h0 = max( h_min, h(k-2) )
95 h1 = max( h_min, h(k-1) )
96 h2 = max( h_min, h(k) )
97 h3 = max( h_min, h(k+1) )
98 endif
99
100 f1 = (h0+h1) * (h2+h3) / (h1+h2)
101 f2 = h2 * u(k-1) + h1 * u(k)
102 f3 = 1.0 / (h0+h1+h2) + 1.0 / (h1+h2+h3)
103 et1 = f1 * f2 * f3
104 et2 = ( h2 * (h2+h3) / ( (h0+h1+h2)*(h0+h1) ) ) * &
105 ((h0+2.0*h1) * u(k-1) - h1 * u(k-2))
106 et3 = ( h1 * (h0+h1) / ( (h1+h2+h3)*(h2+h3) ) ) * &
107 ((2.0*h2+h3) * u(k) - h2 * u(k+1))
108 edge_values(k,1) = (et1 + et2 + et3) / ( h0 + h1 + h2 + h3)
109 edge_values(k-1,2) = edge_values(k,1)
110
111 enddo ! end loop on interior cells
112
113 ! Determine first two edge values
114 h_min = max( this%h_neglect, hminfrac*sum(h(1:4)) )
115 x(1) = 0.0
116 do k = 1,4
117 dx = max(h_min, h(k) )
118 x(k+1) = x(k) + dx
119 do j = 1,4 ; a(k,j) = ( (x(k+1)**j) - (x(k)**j) ) / real(j) ; enddo
120 b(k) = u(k) * dx
121 enddo
122
123 call solve_linear_system( a, b, c, 4 )
124
125 ! Set the edge values of the first cell
126 f = 0.0
127 do k = 1, 4
128 f = f + c(k) * ( x(1)**(k-1) )
129 enddo
130 edge_values(1,1) = f
131 f = 0.0
132 do k = 1, 4
133 f = f + c(k) * ( x(2)**(k-1) )
134 enddo
135 edge_values(1,2) = f
136 edge_values(2,1) = edge_values(1,2)
137
138 ! Determine two edge values of the last cell
139 h_min = max( this%h_neglect, hminfrac*sum(h(n-3:n)) )
140 x(1) = 0.0
141 do k = 1,4
142 dx = max(h_min, h(n-4+k) )
143 x(k+1) = x(k) + dx
144 do j = 1,4 ; a(k,j) = ( (x(k+1)**j) - (x(k)**j) ) / real(j) ; enddo
145 b(k) = u(n-4+k) * dx
146 enddo
147
148 call solve_linear_system( a, b, c, 4 )
149
150 ! Set the last and second to last edge values
151 f = 0.0
152 do k = 1, 4
153 f = f + c(k) * ( x(5)**(k-1) )
154 enddo
155 edge_values(n,2) = f
156 f = 0.0
157 do k = 1, 4
158 f = f + c(k) * ( x(4)**(k-1) )
159 enddo
160 edge_values(n,1) = f
161 edge_values(n-1,2) = edge_values(n,1)
162
163 ! Bound edge values
164 call bound_edge_values( n, h, u, edge_values, this%h_neglect, answer_date=20180101 )
165
166 ! Make discontinuous edge values monotonic
167 call check_discontinuous_edge_values( n, u, edge_values )
168
169 ! Loop on interior cells to apply the standard
170 ! PPM limiter (Colella & Woodward, JCP 84)
171 do k = 2,n-1
172
173 ! Get cell averages
174 u_l = u(k-1)
175 u_c = u(k)
176 u_r = u(k+1)
177
178 edge_l = edge_values(k,1)
179 edge_r = edge_values(k,2)
180
181 if ( (u_r - u_c)*(u_c - u_l) <= 0.0) then
182 ! Flatten extremum
183 edge_l = u_c
184 edge_r = u_c
185 else
186 expr1 = 3.0 * (edge_r - edge_l) * ( (u_c - edge_l) + (u_c - edge_r))
187 expr2 = (edge_r - edge_l) * (edge_r - edge_l)
188 if ( expr1 > expr2 ) then
189 ! Place extremum at right edge of cell by adjusting left edge value
190 edge_l = u_c + 2.0 * ( u_c - edge_r )
191 edge_l = max( min( edge_l, max(u_l, u_c) ), min(u_l, u_c) ) ! In case of round off
192 elseif ( expr1 < -expr2 ) then
193 ! Place extremum at left edge of cell by adjusting right edge value
194 edge_r = u_c + 2.0 * ( u_c - edge_l )
195 edge_r = max( min( edge_r, max(u_r, u_c) ), min(u_r, u_c) ) ! In case of round off
196 endif
197 endif
198 ! This checks that the difference in edge values is representable
199 ! and avoids overshoot problems due to round off.
200 !### The 1.e-60 needs to have units of [A], so this dimensionally inconsistent.
201 if ( abs( edge_r - edge_l )<max(1.e-60,epsilon(u_c)*abs(u_c)) ) then
202 edge_l = u_c
203 edge_r = u_c
204 endif
205
206 edge_values(k,1) = edge_l
207 edge_values(k,2) = edge_r
208
209 enddo ! end loop on interior cells
210
211 ! PCM within boundary cells
212 edge_values(1,:) = u(1)
213 edge_values(n,:) = u(n)
214
215 ! Store reconstruction
216 do k = 1, n
217 this%u_mean(k) = u(k)
218 this%ul(k) = edge_values(k,1)
219 this%ur(k) = edge_values(k,2)
220 enddo
221
222end subroutine reconstruct
223
224!> Runs PPM_H4_2018 reconstruction unit tests and returns True for any fails, False otherwise
225logical function unit_tests(this, verbose, stdout, stderr)
226 class(ppm_h4_2018), intent(inout) :: this !< This reconstruction
227 logical, intent(in) :: verbose !< True, if verbose
228 integer, intent(in) :: stdout !< I/O channel for stdout
229 integer, intent(in) :: stderr !< I/O channel for stderr
230 ! Local variables
231 real, allocatable :: ul(:), ur(:), um(:) ! test values [A]
232 real, allocatable :: ull(:), urr(:) ! test values [A]
233 type(testing) :: test ! convenience functions
234 integer :: k
235
236 call test%set( stdout=stdout ) ! Sets the stdout channel in test
237 call test%set( stderr=stderr ) ! Sets the stderr channel in test
238 call test%set( verbose=verbose ) ! Sets the verbosity flag in test
239
240 if (verbose) write(stdout,'(a)') 'PPM_H4_2018:unit_tests testing with linear fn'
241
242 call this%init(5)
243 call test%test( this%n /= 5, 'Setting number of levels')
244 allocate( um(5), ul(5), ur(5), ull(5), urr(5) )
245
246 ! Straight line, f(x) = x , or f(K) = 2*K
247 call this%reconstruct( (/2.,2.,2.,2.,2./), (/1.,3.,5.,7.,9./) )
248 call test%real_arr(5, this%u_mean, (/1.,3.,5.,7.,9./), 'Setting cell values')
249 call test%real_arr(5, this%ul, (/1.,2.,4.,6.,9./), 'Left edge values', robits=2)
250 call test%real_arr(5, this%ur, (/1.,4.,6.,8.,9./), 'Right edge values', robits=1)
251 do k = 1, 5
252 um(k) = this%u_mean(k)
253 enddo
254 call test%real_arr(5, um, (/1.,3.,5.,7.,9./), 'Return cell mean')
255
256 do k = 1, 5
257 ul(k) = this%f(k, 0.)
258 um(k) = this%f(k, 0.5)
259 ur(k) = this%f(k, 1.)
260 enddo
261 call test%real_arr(5, ul, this%ul, 'Evaluation on left edge')
262 call test%real_arr(5, um, (/1.,3.,5.,7.,9./), 'Evaluation in center')
263 call test%real_arr(5, ur, this%ur, 'Evaluation on right edge')
264
265 do k = 1, 5
266 ul(k) = this%dfdx(k, 0.)
267 um(k) = this%dfdx(k, 0.5)
268 ur(k) = this%dfdx(k, 1.)
269 enddo
270 call test%real_arr(5, ul, (/0.,2.,2.,2.,0./), 'dfdx on left edge', robits=4)
271 call test%real_arr(5, um, (/0.,2.,2.,2.,0./), 'dfdx in center', robits=2)
272 call test%real_arr(5, ur, (/0.,2.,2.,2.,0./), 'dfdx on right edge', robits=6)
273
274 do k = 1, 5
275 um(k) = this%average(k, 0.5, 0.75) ! Average from x=0.25 to 0.75 in each cell
276 enddo
277 call test%real_arr(5, um, (/1.,3.25,5.25,7.25,9./), 'Return interval average')
278
279 if (verbose) write(stdout,'(a)') 'PPM_H4_2018:unit_tests testing with parabola'
280
281 ! x = 3 i i=0 at origin
282 ! f(x) = x^2 / 3 = 3 i^2
283 ! f[i] = [ ( 3 i )^3 - ( 3 i - 3 )^3 ] i=1,2,3,4,5
284 ! means: 1, 7, 19, 37, 61
285 ! edges: 0, 3, 12, 27, 48, 75
286 call this%reconstruct( (/3.,3.,3.,3.,3./), (/1.,7.,19.,37.,61./) )
287 do k = 1, 5
288 ul(k) = this%f(k, 0.)
289 um(k) = this%f(k, 0.5)
290 ur(k) = this%f(k, 1.)
291 enddo
292 call test%real_arr(5, ul, (/1.,3.,12.,27.,61./), 'Return left edge', robits=2)
293 call test%real_arr(5, ur, (/1.,12.,27.,48.,61./), 'Return right edge', robits=1)
294
295 call this%destroy()
296 deallocate( um, ul, ur, ull, urr )
297
298 unit_tests = test%summarize('PPM_H4_2018:unit_tests')
299
300end function unit_tests
301
302!> \namespace recon1d_ppm_h4_2018
303!!
304
305end module recon1d_ppm_h4_2018