Recon1d_PPM_hybgen.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 following Colella and Woodward, 1984
6!!
7!! This implementation of PPM follows Colella and Woodward, 1984 \cite colella1984, with
8!! cells resorting to PCM for extrema including first and last cells in column. The algorithm was
9!! first ported from Hycom as hybgen_ppm_coefs() in the mom_hybgen_remap module. This module is
10!! a refactor to facilitate more complete testing and evaluation.
11!!
12!! The mom_hybgen_remap.hybgen_ppm_coefs() function (reached with "PPM_HYGEN"),
13!! regrid_edge_values.edge_values_explicit_h4cw() function followed by ppm_functions.ppm_reconstruction()
14!! (reached with "PPM_CW"), are equivalent. Similarly recon1d_ppm_hybgen (this implementation) is equivalent also.
16
17use recon1d_type, only : testing
18use recon1d_ppm_cw, only : ppm_cw
19
20implicit none ; private
21
22public ppm_hybgen, testing
23
24!> PPM reconstruction following White and Adcroft, 2008
25!!
26!! Implemented by extending recon1d_ppm_cwk.
27!!
28!! The source for the methods ultimately used by this class are:
29!! - init() -> recon1d_ppm_cw.init()
30!! - reconstruct() *locally defined
31!! - average() -> recon1d_ppm_cw.average()
32!! - f() -> recon1d_ppm_cw.f()
33!! - dfdx() -> recon1d_ppm_cw.dfdx()
34!! - check_reconstruction() *locally defined
35!! - unit_tests() -> recon1d_ppm_cw.unit_tests()
36!! - destroy() -> recon1d_ppm_cw.destroy()
37!! - remap_to_sub_grid() -> recon1d_type.remap_to_sub_grid()
38!! - init_parent() -> init()
39!! - reconstruct_parent() -> reconstruct()
40type, extends (ppm_cw) :: ppm_hybgen
41
42contains
43 !> Implementation of the PPM_hybgen reconstruction
44 procedure :: reconstruct => reconstruct
45 !> Implementation of check reconstruction for the PPM_hybgen reconstruction
47 !> Implementation of unit tests for the PPM_hybgen reconstruction
48 procedure :: unit_tests => unit_tests
49
50end type ppm_hybgen
51
52contains
53
54!> Calculate a 1D PPM_hybgen reconstructions based on h(:) and u(:)
55subroutine reconstruct(this, h, u)
56 class(ppm_hybgen), intent(inout) :: this !< This reconstruction
57 real, intent(in) :: h(*) !< Grid spacing (thickness) [typically H]
58 real, intent(in) :: u(*) !< Cell mean values [A]
59 ! Local variables
60 real :: h0, h1, h2, h3 ! Cell thickness h(k-2), h(k-1), h(k), h(k+1) in K loop [H]
61 real :: h01_h112, h23_h122 ! Approximately 2/3 [nondim]
62 real :: h112, h122 ! Approximately 3 h [H]
63 real :: ddh ! Approximately 0 [nondim]
64 real :: I_h12, I_h01, I_h0123 ! Reciprocals of d12 and sum(h) [H-1]
65 real :: dul, dur ! Left and right cell PLM slopes [A]
66 real :: u0, u1, u2 ! Far left, left, and right cell values [A]
67 real :: edge ! Edge value between cell k-1 and k [A]
68 real :: u_min, u_max ! Minimum and maximum value across edge [A]
69 real :: a6 ! Colella and Woodward curvature [A]
70 real :: du, duc ! Difference between edges across cell [A]
71 real :: slp(this%n) ! PLM slope [A]
72 real :: sigma_l, sigma_c, sigma_r ! Left, central and right slope estimates as
73 ! differences across the cell [A]
74 real :: slope_x_h ! retained PLM slope times half grid step [A]
75 real :: edge_l, edge_r ! Edge values (left and right) [A]
76 real :: expr1, expr2 ! Temporary expressions [A2]
77 real :: u0_avg ! avg value at given edge [A]
78 integer :: k, n, km1, kp1
79
80 n = this%n
81
82 ! First populate the PLM reconstructions
83 slp(1) = 0.
84 do k = 2, n-1
85 h0 = max( this%h_neglect, h(k-1) )
86 h1 = max( this%h_neglect, h(k) )
87 h2 = max( this%h_neglect, h(k+1) )
88 dul = u(k) - u(k-1)
89 dur = u(k+1) - u(k)
90 h112 = ( 2.0 * h0 + h1 )
91 h122 = ( h1 + 2.0 * h2 )
92 i_h01 = 1. / ( h0 + h1 )
93 i_h12 = 1. / ( h1 + h2 )
94 h01_h112 = ( 2.0 * h0 + h1 ) / ( h0 + h1 ) ! When uniform -> 3/2
95 h23_h122 = ( 2.0 * h2 + h1 ) / ( h2 + h1 ) ! When uniform -> 3/2
96 if ( dul * dur > 0.) then
97 du = ( h1 / ( h1 + ( h0 + h2 ) ) ) * ( h112 * dur * i_h12 + h122 * dul * i_h01 )
98 slp(k) = sign( min( abs(2.0 * dul), abs(du), abs(2.0 * dur) ), du)
99 else
100 slp(k) = 0.
101 endif
102 enddo
103 slp(n) = 0.
104
105 this%ul(1) = u(1) ! PCM
106 this%ur(1) = u(1) ! PCM
107 this%ul(2) = u(1) ! PCM
108 do k = 3, n-1 ! K=3 is interface between cells 2 and 3
109 h0 = max( this%h_neglect, h(k-2) )
110 h1 = max( this%h_neglect, h(k-1) )
111 h2 = max( this%h_neglect, h(k) )
112 h3 = max( this%h_neglect, h(k+1) )
113 h01_h112 = ( h0 + h1 ) / ( 2. * h1 + h2 ) ! When uniform -> 2/3
114 h23_h122 = ( h2 + h3 ) / ( h1 + 2. * h2 ) ! When uniform -> 2/3
115 ddh = h01_h112 - h23_h122 ! When uniform -> 0
116 i_h12 = 1.0 / ( h1 + h2 ) ! When uniform -> 1/(2h)
117 i_h0123 = 1.0 / ( ( h0 + h1 ) + ( h2 + h3 ) ) ! When uniform -> 1/(4h)
118 dul = slp(k-1)
119 dur = slp(k)
120 u1 = u(k-1)
121 u2 = u(k)
122 edge = i_h12 * ( h2 * u1 + h1 * u2 ) & ! 1/2 u1 + 1/2 u2
123 + i_h0123 * ( 2.0 * h1 * h2 * i_h12 * ( u2 - u1 ) * ddh & ! 0
124 + ( h2 * dul * h23_h122 - h1 * dur * h01_h112 ) ) ! 1/6 dul - 1/6 dur
125 this%ur(k-1) = edge
126 this%ul(k) = edge
127 enddo
128 this%ur(n-1) = u(n) ! PCM
129 this%ur(n) = u(n) ! PCM
130 this%ul(n) = u(n) ! PCM
131
132 do k = 2, n-1 ! K=2 is interface between cells 1 and 2
133 u0 = u(k-1)
134 u1 = u(k)
135 u2 = u(k+1)
136 a6 = 3.0 * ( ( u1 - this%ul(k) ) + ( u1 - this%ur(k) ) )
137 a6 = 6.0 * u1 - 3.0 * ( this%ul(k) + this%ur(k) )
138 du = this%ur(k) - this%ul(k)
139 if ( ( u2 - u1 ) * ( u1 - u0 ) <= 0.0 ) then ! Large scale extrema
140 this%ul(k) = u1
141 this%ur(k) = u1
142 elseif ( du * a6 > du * du ) then ! Extrema on right
143 edge = 3.0 * u1 - 2.0 * this%ur(k) ! Subject to round off
144 ! u_min = min( u0, u1 )
145 ! u_max = max( u0, u1 )
146 ! edge = max( min( edge, u_max), u_min )
147 this%ul(k) = edge
148 elseif ( du * a6 < - du * du ) then ! Extrema on left
149 edge = 3.0 * u1 - 2.0 * this%ul(k) ! Subject to round off
150 ! u_min = min( u1, u2 )
151 ! u_max = max( u1, u2 )
152 ! edge = max( min( edge, u_max), u_min )
153 this%ur(k) = edge
154 endif
155 enddo
156
157 ! ### Note that the PPM_HYBGEM option calculated the CW PPM coefficients and then
158 ! invoked the OM4-era limiters afterwards, effectively doing the limiters twice.
159 ! This second pass does change answers!
160
161 ! Loop on cells to bound edge value
162 do k = 1, n
163
164 ! For the sake of bounding boundary edge values, the left neighbor of the left boundary cell
165 ! is assumed to be the same as the left boundary cell and the right neighbor of the right
166 ! boundary cell is assumed to be the same as the right boundary cell. This effectively makes
167 ! boundary cells look like extrema.
168 km1 = max(1,k-1) ; kp1 = min(k+1,n)
169
170 slope_x_h = 0.0
171 sigma_l = ( u(k) - u(km1) )
172 if ( (h(km1) + h(kp1)) + 2.0*h(k) > 0. ) then
173 sigma_c = ( u(kp1) - u(km1) ) * ( h(k) / ((h(km1) + h(kp1)) + 2.0*h(k)) )
174 else
175 sigma_c = 0.
176 endif
177 sigma_r = ( u(kp1) - u(k) )
178
179 ! The limiter is used in the local coordinate system to each cell, so for convenience store
180 ! the slope times a half grid spacing. (See White and Adcroft JCP 2008 Eqs 19 and 20)
181 if ( (sigma_l * sigma_r) > 0.0 ) &
182 slope_x_h = sign( min(abs(sigma_l),abs(sigma_c),abs(sigma_r)), sigma_c )
183
184 ! Limit the edge values
185 if ( (u(km1)-this%ul(k)) * (this%ul(k)-u(k)) < 0.0 ) then
186 this%ul(k) = u(k) - sign( min( abs(slope_x_h), abs(this%ul(k)-u(k)) ), slope_x_h )
187 endif
188
189 if ( (u(kp1)-this%ur(k)) * (this%ur(k)-u(k)) < 0.0 ) then
190 this%ur(k) = u(k) + sign( min( abs(slope_x_h), abs(this%ur(k)-u(k)) ), slope_x_h )
191 endif
192
193 ! Finally bound by neighboring cell means in case of roundoff
194 this%ul(k) = max( min( this%ul(k), max(u(km1), u(k)) ), min(u(km1), u(k)) )
195 this%ur(k) = max( min( this%ur(k), max(u(kp1), u(k)) ), min(u(kp1), u(k)) )
196
197 enddo ! loop on interior edges
198
199 do k = 1, n-1
200 if ( (this%ul(k+1) - this%ur(k)) * (u(k+1) - u(k)) < 0.0 ) then
201 u0_avg = 0.5 * ( this%ur(k) + this%ul(k+1) )
202 u0_avg = max( min( u0_avg, max(u(k), u(k+1)) ), min(u(k), u(k+1)) )
203 this%ur(k) = u0_avg
204 this%ul(k+1) = u0_avg
205 endif
206 enddo ! end loop on interior edges
207
208 ! Loop on interior cells to apply the standard
209 ! PPM limiter (Colella & Woodward, JCP 84)
210 do k = 2, n-1
211
212 ! Get cell averages
213 u0 = u(k-1)
214 u1 = u(k)
215 u2 = u(k+1)
216
217 edge_l = this%ul(k)
218 edge_r = this%ur(k)
219
220 if ( (u2 - u1)*(u1 - u0) <= 0.0) then
221 ! Flatten extremum
222 edge_l = u1
223 edge_r = u1
224 else
225 expr1 = 3.0 * (edge_r - edge_l) * ( (u1 - edge_l) + (u1 - edge_r))
226 expr2 = (edge_r - edge_l) * (edge_r - edge_l)
227 if ( expr1 > expr2 ) then
228 ! Place extremum at right edge of cell by adjusting left edge value
229 edge_l = u1 + 2.0 * ( u1 - edge_r )
230 edge_l = max( min( edge_l, max(u0, u1) ), min(u0, u1) ) ! In case of round off
231 elseif ( expr1 < -expr2 ) then
232 ! Place extremum at left edge of cell by adjusting right edge value
233 edge_r = u1 + 2.0 * ( u1 - edge_l )
234 edge_r = max( min( edge_r, max(u2, u1) ), min(u2, u1) ) ! In case of round off
235 endif
236 endif
237 ! This checks that the difference in edge values is representable
238 ! and avoids overshoot problems due to round off.
239 !### The 1.e-60 needs to have units of [A], so this dimensionally inconsistent.
240 if ( abs( edge_r - edge_l )<max(1.e-60,epsilon(u1)*abs(u1)) ) then
241 edge_l = u1
242 edge_r = u1
243 endif
244
245 this%ul(k) = edge_l
246 this%ur(k) = edge_r
247
248 enddo ! end loop on interior cells
249
250
251 ! After the limiter, are ur and ul bounded???? -AJA
252
253 ! Store mean
254 do k = 1, n
255 this%u_mean(k) = u(k)
256 enddo
257
258end subroutine reconstruct
259
260!> Checks the PPM_hybgen reconstruction for consistency
261logical function check_reconstruction(this, h, u)
262 class(ppm_hybgen), intent(in) :: this !< This reconstruction
263 real, intent(in) :: h(*) !< Grid spacing (thickness) [typically H]
264 real, intent(in) :: u(*) !< Cell mean values [A]
265 ! Local variables
266 integer :: k
267
268 check_reconstruction = .false.
269
270 ! Simply checks the internal copy of "u" is exactly equal to "u"
271 do k = 1, this%n
272 if ( abs( this%u_mean(k) - u(k) ) > 0. ) check_reconstruction = .true.
273 enddo
274
275 ! If (u - ul) has the opposite sign from (ur - u), then this cell has an interior extremum
276 do k = 1, this%n
277 if ( ( this%u_mean(k) - this%ul(k) ) * ( this%ur(k) - this%u_mean(k) ) < 0. ) check_reconstruction = .true.
278 enddo
279
280 ! The following consistency checks would fail for this implementation of PPM CW,
281 ! due to round off in the final limiter violating the monotonicity of edge values,
282 ! but actually passes due to the second pass of the limiters with explicit bounding.
283 ! i.e. This implementation cheats!
284
285 ! Check bounding of right edges, w.r.t. the cell means
286 do k = 1, this%n-1
287 if ( ( this%ur(k) - this%u_mean(k) ) * ( this%u_mean(k+1) - this%ur(k) ) < 0. ) check_reconstruction = .true.
288 enddo
289
290 ! Check bounding of left edges, w.r.t. the cell means
291 do k = 2, this%n
292 if ( ( this%u_mean(k) - this%ul(k) ) * ( this%ul(k) - this%u_mean(k-1) ) < 0. ) check_reconstruction = .true.
293 enddo
294
295 ! Check bounding of right edges, w.r.t. this cell mean and the next cell left edge
296 do k = 1, this%n-1
297 if ( ( this%ur(k) - this%u_mean(k) ) * ( this%ul(k+1) - this%ur(k) ) < 0. ) check_reconstruction = .true.
298 enddo
299
300 ! Check bounding of left edges, w.r.t. this cell mean and the previous cell right edge
301 do k = 2, this%n
302 if ( ( this%u_mean(k) - this%ul(k) ) * ( this%ul(k) - this%ur(k-1) ) < 0. ) check_reconstruction = .true.
303 enddo
304
305end function check_reconstruction
306
307!> Runs PPM_hybgen reconstruction unit tests and returns True for any fails, False otherwise
308logical function unit_tests(this, verbose, stdout, stderr)
309 class(ppm_hybgen), intent(inout) :: this !< This reconstruction
310 logical, intent(in) :: verbose !< True, if verbose
311 integer, intent(in) :: stdout !< I/O channel for stdout
312 integer, intent(in) :: stderr !< I/O channel for stderr
313 ! Local variables
314 real, allocatable :: ul(:), ur(:), um(:) ! test values [A]
315 real, allocatable :: ull(:), urr(:) ! test values [A]
316 type(testing) :: test ! convenience functions
317 integer :: k
318
319 call test%set( stdout=stdout ) ! Sets the stdout channel in test
320 call test%set( stderr=stderr ) ! Sets the stderr channel in test
321 call test%set( verbose=verbose ) ! Sets the verbosity flag in test
322
323 if (verbose) write(stdout,'(a)') 'PPM_hybgen:unit_tests testing with linear fn'
324
325 call this%init(5)
326 call test%test( this%n /= 5, 'Setting number of levels')
327 allocate( um(5), ul(5), ur(5), ull(5), urr(5) )
328
329 ! Straight line, f(x) = x , or f(K) = 2*K
330 call this%reconstruct( (/2.,2.,2.,2.,2./), (/1.,4.,7.,10.,13./) )
331 call test%real_arr(5, this%u_mean, (/1.,4.,7.,10.,13./), 'Setting cell values')
332 ! Without PLM extrapolation we get l(2)=2 and r(4)=12 due to PLM=0 in boundary cells. -AJA
333 call test%real_arr(5, this%ul, (/1.,1.,5.5,8.5,13./), 'Left edge values')
334 call test%real_arr(5, this%ur, (/1.,5.5,8.5,13.,13./), 'Right edge values')
335
336 do k = 1, 5
337 ul(k) = this%f(k, 0.)
338 um(k) = this%f(k, 0.5)
339 ur(k) = this%f(k, 1.)
340 enddo
341 call test%real_arr(5, ul, this%ul, 'Evaluation on left edge')
342 call test%real_arr(5, um, (/1.,4.375,7.,9.625,13./), 'Evaluation in center')
343 call test%real_arr(5, ur, this%ur, 'Evaluation on right edge')
344
345 do k = 1, 5
346 ul(k) = this%dfdx(k, 0.)
347 um(k) = this%dfdx(k, 0.5)
348 ur(k) = this%dfdx(k, 1.)
349 enddo
350 ! Most of these values are affected by the PLM boundary cells
351 call test%real_arr(5, ul, (/0.,0.,3.,9.,0./), 'dfdx on left edge')
352 call test%real_arr(5, um, (/0.,4.5,3.,4.5,0./), 'dfdx in center')
353 call test%real_arr(5, ur, (/0.,9.,3.,0.,0./), 'dfdx on right edge')
354
355 do k = 1, 5
356 um(k) = this%average(k, 0.5, 0.75) ! Average from x=0.25 to 0.75 in each cell
357 enddo
358 ! Most of these values are affected by the PLM boundary cells
359 call test%real_arr(5, um, (/1.,4.84375,7.375,10.28125,13./), 'Return interval average')
360
361 if (verbose) write(stdout,'(a)') 'PPM_hybgen:unit_tests testing with parabola'
362
363 ! x = 2 i i=0 at origin
364 ! f(x) = 3/4 x^2 = (2 i)^2
365 ! f[i] = 3/4 ( 2 i - 1 )^2 on centers
366 ! f[I] = 3/4 ( 2 I )^2 on edges
367 ! f[i] = 1/8 [ x^3 ] for means
368 ! edges: 0, 1, 12, 27, 48, 75
369 ! means: 1, 7, 19, 37, 61
370 ! cengters: 0.75, 6.75, 18.75, 36.75, 60.75
371 call this%reconstruct( (/2.,2.,2.,2.,2./), (/1.,7.,19.,37.,61./) )
372 do k = 1, 5
373 ul(k) = this%f(k, 0.)
374 um(k) = this%f(k, 0.5)
375 ur(k) = this%f(k, 1.)
376 enddo
377 call test%real_arr(5, ul, (/1.,1.,12.,27.,61./), 'Return left edge')
378 call test%real_arr(5, um, (/1.,7.25,18.75,34.5,61./), 'Return center')
379 call test%real_arr(5, ur, (/1.,12.,27.,57.,61./), 'Return right edge')
380
381 ! x = 3 i i=0 at origin
382 ! f(x) = x^2 / 3 = 3 i^2
383 ! f[i] = [ ( 3 i )^3 - ( 3 i - 3 )^3 ] i=1,2,3,4,5
384 ! means: 1, 7, 19, 37, 61
385 ! edges: 0, 3, 12, 27, 48, 75
386 call this%reconstruct( (/3.,3.,3.,3.,3./), (/1.,7.,19.,37.,61./) )
387 do k = 1, 5
388 ul(k) = this%f(k, 0.)
389 um(k) = this%f(k, 0.5)
390 ur(k) = this%f(k, 1.)
391 enddo
392 call test%real_arr(5, ul, (/1.,1.,12.,27.,61./), 'Return left edge')
393 call test%real_arr(5, ur, (/1.,12.,27.,57.,61./), 'Return right edge')
394
395 call this%destroy()
396 deallocate( um, ul, ur, ull, urr )
397
398 unit_tests = test%summarize('PPM_hybgen:unit_tests')
399
400end function unit_tests
401
402!> \namespace recon1d_ppm_hybgen
403!!
404
405end module recon1d_ppm_hybgen