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() *locally defined but calls 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 the PPM_hybgen average over an interval [A]
46 procedure :: average => average
47 !> Implementation of check reconstruction for the PPM_hybgen reconstruction
49 !> Implementation of unit tests for the PPM_hybgen reconstruction
50 procedure :: unit_tests => unit_tests
51
52end type ppm_hybgen
53
54contains
55
56!> Calculate a 1D PPM_hybgen reconstruction based on h(:) and u(:)
57!!
58!! First pass: hybgen_ppm_coefs() computes initial edge estimates with CW monotonicity.
59!! Second pass: applies OM4-era bound_edge_values() and check_discontinuous_edge_values(),
60!! then the standard CW PPM limiter (post-2018 expressions, answer_date=99991231).
61!! This reproduces bit-for-bit the behavior of the old-style PPM_HYBGEN scheme.
62subroutine reconstruct(this, h, u)
63 class(ppm_hybgen), intent(inout) :: this !< This reconstruction
64 real, intent(in) :: h(*) !< Grid spacing (thickness) [typically H]
65 real, intent(in) :: u(*) !< Cell mean values [A]
66 ! Local variables
67 integer :: k, n
68 real :: ppoly_e(this%n, 2) ! PPM edge values [A]
69 real :: u_l, u_c, u_r ! Left, center, right cell averages [A]
70 real :: edge_l, edge_r ! Left and right edge values [A]
71 real :: expr1, expr2 ! Temporary expressions [A2]
72
73 n = this%n
74
75 ! First pass: compute initial edge estimates using the hybgen algorithm with CW monotonicity
76 call hybgen_ppm_coefs(u, h, ppoly_e, n, this%h_neglect)
77
78 ! Second pass: apply OM4-era PPM limiters (post-2018 answers via answer_date=99991231)
79 call bound_edge_values(n, h, u, ppoly_e, this%h_neglect, answer_date=99991231)
80 call check_discontinuous_edge_values(n, u, ppoly_e)
81
82 ! Apply the standard CW PPM limiter (Colella & Woodward, JCP 84) on interior cells
83 do k = 2, n-1
84 u_l = u(k-1) ; u_c = u(k) ; u_r = u(k+1)
85 edge_l = ppoly_e(k,1) ; edge_r = ppoly_e(k,2)
86 if ( (u_r - u_c)*(u_c - u_l) <= 0.0 ) then
87 edge_l = u_c ; edge_r = u_c
88 else
89 expr1 = 3.0 * (edge_r - edge_l) * ( (u_c - edge_l) + (u_c - edge_r) )
90 expr2 = (edge_r - edge_l) * (edge_r - edge_l)
91 if ( expr1 > expr2 ) then
92 edge_l = u_c + 2.0 * ( u_c - edge_r )
93 edge_l = max( min( edge_l, max(u_l, u_c) ), min(u_l, u_c) )
94 elseif ( expr1 < -expr2 ) then
95 edge_r = u_c + 2.0 * ( u_c - edge_l )
96 edge_r = max( min( edge_r, max(u_r, u_c) ), min(u_r, u_c) )
97 endif
98 endif
99 !### The 1.e-60 needs to have units of [A], so this is dimensionally inconsistent.
100 if ( abs( edge_r - edge_l ) < max(1.e-60, epsilon(u_c)*abs(u_c)) ) then
101 edge_l = u_c ; edge_r = u_c
102 endif
103 ppoly_e(k,1) = edge_l ; ppoly_e(k,2) = edge_r
104 enddo
105 ! Boundary cells are PCM
106 ppoly_e(1,:) = u(1) ; ppoly_e(n,:) = u(n)
107
108 do k = 1, n
109 this%ul(k) = ppoly_e(k, 1)
110 this%ur(k) = ppoly_e(k, 2)
111 this%u_mean(k) = u(k)
112 enddo
113
114end subroutine reconstruct
115
116!> Average between xa and xb for cell k of a PPM_hybgen reconstruction [A]
117!!
118!! Calls the parent PPM_CW average and then clamps the result to [min(ul,ur), max(ul,ur)].
119!! This replicates the force_bounds_in_subcell behavior of the equivalent old-style PPM_HYBGEN
120!! scheme.
121real function average(this, k, xa, xb)
122 class(ppm_hybgen), intent(in) :: this !< This reconstruction
123 integer, intent(in) :: k !< Cell number
124 real, intent(in) :: xa !< Start of averaging interval on element (0 to 1)
125 real, intent(in) :: xb !< End of averaging interval on element (0 to 1)
126 real :: u_lo, u_hi ! Bounds on the sub-cell average given by the edge values [A]
127
128 average = this%PPM_CW%average(k, xa, xb)
129 u_lo = min(this%ul(k), this%ur(k))
130 u_hi = max(this%ul(k), this%ur(k))
131 average = max(u_lo, min(u_hi, average))
132
133end function average
134
135!> Checks the PPM_hybgen reconstruction for consistency
136logical function check_reconstruction(this, h, u)
137 class(ppm_hybgen), intent(in) :: this !< This reconstruction
138 real, intent(in) :: h(*) !< Grid spacing (thickness) [typically H]
139 real, intent(in) :: u(*) !< Cell mean values [A]
140 ! Local variables
141 integer :: k
142
143 check_reconstruction = .false.
144
145 ! Simply checks the internal copy of "u" is exactly equal to "u"
146 do k = 1, this%n
147 if ( abs( this%u_mean(k) - u(k) ) > 0. ) check_reconstruction = .true.
148 enddo
149
150 ! If (u - ul) has the opposite sign from (ur - u), then this cell has an interior extremum
151 do k = 1, this%n
152 if ( ( this%u_mean(k) - this%ul(k) ) * ( this%ur(k) - this%u_mean(k) ) < 0. ) check_reconstruction = .true.
153 enddo
154
155 ! The following consistency checks would fail for this implementation of PPM CW,
156 ! due to round off in the final limiter violating the monotonicity of edge values,
157 ! but actually passes due to the second pass of the limiters with explicit bounding.
158 ! i.e. This implementation cheats!
159
160 ! Check bounding of right edges, w.r.t. the cell means
161 do k = 1, this%n-1
162 if ( ( this%ur(k) - this%u_mean(k) ) * ( this%u_mean(k+1) - this%ur(k) ) < 0. ) check_reconstruction = .true.
163 enddo
164
165 ! Check bounding of left edges, w.r.t. the cell means
166 do k = 2, this%n
167 if ( ( this%u_mean(k) - this%ul(k) ) * ( this%ul(k) - this%u_mean(k-1) ) < 0. ) check_reconstruction = .true.
168 enddo
169
170 ! Check bounding of right edges, w.r.t. this cell mean and the next cell left edge
171 do k = 1, this%n-1
172 if ( ( this%ur(k) - this%u_mean(k) ) * ( this%ul(k+1) - this%ur(k) ) < 0. ) check_reconstruction = .true.
173 enddo
174
175 ! Check bounding of left edges, w.r.t. this cell mean and the previous cell right edge
176 do k = 2, this%n
177 if ( ( this%u_mean(k) - this%ul(k) ) * ( this%ul(k) - this%ur(k-1) ) < 0. ) check_reconstruction = .true.
178 enddo
179
180end function check_reconstruction
181
182!> Runs PPM_hybgen reconstruction unit tests and returns True for any fails, False otherwise
183logical function unit_tests(this, verbose, stdout, stderr)
184 class(ppm_hybgen), intent(inout) :: this !< This reconstruction
185 logical, intent(in) :: verbose !< True, if verbose
186 integer, intent(in) :: stdout !< I/O channel for stdout
187 integer, intent(in) :: stderr !< I/O channel for stderr
188 ! Local variables
189 real, allocatable :: ul(:), ur(:), um(:) ! test values [A]
190 real, allocatable :: ull(:), urr(:) ! test values [A]
191 type(testing) :: test ! convenience functions
192 integer :: k
193
194 call test%set( stdout=stdout ) ! Sets the stdout channel in test
195 call test%set( stderr=stderr ) ! Sets the stderr channel in test
196 call test%set( verbose=verbose ) ! Sets the verbosity flag in test
197
198 if (verbose) write(stdout,'(a)') 'PPM_hybgen:unit_tests testing with linear fn'
199
200 call this%init(5)
201 call test%test( this%n /= 5, 'Setting number of levels')
202 allocate( um(5), ul(5), ur(5), ull(5), urr(5) )
203
204 ! Straight line, f(x) = x , or f(K) = 2*K
205 call this%reconstruct( (/2.,2.,2.,2.,2./), (/1.,4.,7.,10.,13./) )
206 call test%real_arr(5, this%u_mean, (/1.,4.,7.,10.,13./), 'Setting cell values')
207 ! Without PLM extrapolation we get l(2)=2 and r(4)=12 due to PLM=0 in boundary cells. -AJA
208 call test%real_arr(5, this%ul, (/1.,1.,5.5,8.5,13./), 'Left edge values')
209 call test%real_arr(5, this%ur, (/1.,5.5,8.5,13.,13./), 'Right edge values')
210
211 do k = 1, 5
212 ul(k) = this%f(k, 0.)
213 um(k) = this%f(k, 0.5)
214 ur(k) = this%f(k, 1.)
215 enddo
216 call test%real_arr(5, ul, this%ul, 'Evaluation on left edge')
217 call test%real_arr(5, um, (/1.,4.375,7.,9.625,13./), 'Evaluation in center')
218 call test%real_arr(5, ur, this%ur, 'Evaluation on right edge')
219
220 do k = 1, 5
221 ul(k) = this%dfdx(k, 0.)
222 um(k) = this%dfdx(k, 0.5)
223 ur(k) = this%dfdx(k, 1.)
224 enddo
225 ! Most of these values are affected by the PLM boundary cells
226 call test%real_arr(5, ul, (/0.,0.,3.,9.,0./), 'dfdx on left edge')
227 call test%real_arr(5, um, (/0.,4.5,3.,4.5,0./), 'dfdx in center')
228 call test%real_arr(5, ur, (/0.,9.,3.,0.,0./), 'dfdx on right edge')
229
230 do k = 1, 5
231 um(k) = this%average(k, 0.5, 0.75) ! Average from x=0.25 to 0.75 in each cell
232 enddo
233 ! Most of these values are affected by the PLM boundary cells
234 call test%real_arr(5, um, (/1.,4.84375,7.375,10.28125,13./), 'Return interval average')
235
236 if (verbose) write(stdout,'(a)') 'PPM_hybgen:unit_tests testing with parabola'
237
238 ! x = 2 i i=0 at origin
239 ! f(x) = 3/4 x^2 = (2 i)^2
240 ! f[i] = 3/4 ( 2 i - 1 )^2 on centers
241 ! f[I] = 3/4 ( 2 I )^2 on edges
242 ! f[i] = 1/8 [ x^3 ] for means
243 ! edges: 0, 1, 12, 27, 48, 75
244 ! means: 1, 7, 19, 37, 61
245 ! cengters: 0.75, 6.75, 18.75, 36.75, 60.75
246 call this%reconstruct( (/2.,2.,2.,2.,2./), (/1.,7.,19.,37.,61./) )
247 do k = 1, 5
248 ul(k) = this%f(k, 0.)
249 um(k) = this%f(k, 0.5)
250 ur(k) = this%f(k, 1.)
251 enddo
252 call test%real_arr(5, ul, (/1.,1.,12.,27.,61./), 'Return left edge')
253 call test%real_arr(5, um, (/1.,7.25,18.75,34.5,61./), 'Return center')
254 call test%real_arr(5, ur, (/1.,12.,27.,57.,61./), 'Return right edge')
255
256 ! x = 3 i i=0 at origin
257 ! f(x) = x^2 / 3 = 3 i^2
258 ! f[i] = [ ( 3 i )^3 - ( 3 i - 3 )^3 ] i=1,2,3,4,5
259 ! means: 1, 7, 19, 37, 61
260 ! edges: 0, 3, 12, 27, 48, 75
261 call this%reconstruct( (/3.,3.,3.,3.,3./), (/1.,7.,19.,37.,61./) )
262 do k = 1, 5
263 ul(k) = this%f(k, 0.)
264 um(k) = this%f(k, 0.5)
265 ur(k) = this%f(k, 1.)
266 enddo
267 call test%real_arr(5, ul, (/1.,1.,12.,27.,61./), 'Return left edge')
268 call test%real_arr(5, ur, (/1.,12.,27.,57.,61./), 'Return right edge')
269
270 call this%destroy()
271 deallocate( um, ul, ur, ull, urr )
272
273 unit_tests = test%summarize('PPM_hybgen:unit_tests')
274
275end function unit_tests
276
277!> \namespace recon1d_ppm_hybgen
278!!
279
280! ============================================================================
281! Private subroutines copied from phased-out modules to avoid dependencies.
282! These reproduce bit-for-bit the results of the original functions they replace.
283! ============================================================================
284
285!> Set up edge values for PPM reconstruction using the hybgen (HYCOM) algorithm.
286!!
287!! Copied from MOM_hybgen_remap.hybgen_ppm_coefs().
288!! Original code by Tim Campbell (MSU, 2002) and Alan Wallcraft (NRL, 2007).
289subroutine hybgen_ppm_coefs(s, h_src, edges, nk, thin, PCM_lay)
290 integer, intent(in) :: nk !< The number of input layers
291 real, intent(in) :: s(nk) !< The input scalar fields [A]
292 real, intent(in) :: h_src(nk) !< The input grid layer thicknesses [H ~> m or kg m-2]
293 real, intent(out) :: edges(nk,2) !< The PPM interpolation edge values [A]
294 real, intent(in) :: thin !< A negligible layer thickness [H ~> m or kg m-2]
295 logical, optional, intent(in) :: PCM_lay(nk) !< If true for a layer, use PCM remapping
296
297 real :: dp(nk) ! Input grid layer thicknesses, but with a minimum thickness given by thin [H ~> m or kg m-2]
298 logical :: PCM_layer(nk) ! True for layers that should use PCM remapping
299 real :: da ! Difference between the unlimited scalar edge value estimates [A]
300 real :: a6 ! Scalar field differences that are proportional to the curvature [A]
301 real :: slk, srk ! Differences between adjacent cell averages of scalars [A]
302 real :: sck ! Scalar differences across a cell [A]
303 real :: as(nk) ! Scalar field difference across each cell [A]
304 real :: al(nk), ar(nk) ! Scalar field at the left and right edges of a cell [A]
305 real :: h112(nk+1), h122(nk+1) ! Combinations of thicknesses [H ~> m or kg m-2]
306 real :: I_h12(nk+1) ! Inverses of combinations of thicknesses [H-1 ~> m-1 or m2 kg-1]
307 real :: h2_h123(nk) ! A ratio of a layer thickness to the sum of 3 adjacent thicknesses [nondim]
308 real :: I_h0123(nk) ! Inverse of the sum of 4 adjacent thicknesses [H-1 ~> m-1 or m2 kg-1]
309 real :: h01_h112(nk+1) ! A ratio of sums of adjacent thicknesses [nondim]
310 real :: h23_h122(nk+1) ! A ratio of sums of adjacent thicknesses [nondim]
311 integer :: k
312
313 do k=1,nk ; dp(k) = max(h_src(k), thin) ; enddo
314
315 if (present(pcm_lay)) then
316 do k=1,nk ; pcm_layer(k) = (pcm_lay(k) .or. dp(k) <= thin) ; enddo
317 else
318 do k=1,nk ; pcm_layer(k) = (dp(k) <= thin) ; enddo
319 endif
320
321 do k=2,nk
322 h112(k) = 2.*dp(k-1) + dp(k)
323 h122(k) = dp(k-1) + 2.*dp(k)
324 i_h12(k) = 1.0 / (dp(k-1) + dp(k))
325 enddo
326 do k=2,nk-1
327 h2_h123(k) = dp(k) / (dp(k) + (dp(k-1)+dp(k+1)))
328 enddo
329 do k=3,nk-1
330 i_h0123(k) = 1.0 / ((dp(k-2) + dp(k-1)) + (dp(k) + dp(k+1)))
331 h01_h112(k) = (dp(k-2) + dp(k-1)) / (2.0*dp(k-1) + dp(k))
332 h23_h122(k) = (dp(k) + dp(k+1)) / (dp(k-1) + 2.0*dp(k))
333 enddo
334
335 as(1) = 0.
336 do k=2,nk-1
337 if (pcm_layer(k)) then
338 as(k) = 0.0
339 else
340 slk = s(k)-s(k-1)
341 srk = s(k+1)-s(k)
342 if (slk*srk > 0.) then
343 sck = h2_h123(k)*( h112(k)*srk*i_h12(k+1) + h122(k+1)*slk*i_h12(k) )
344 as(k) = sign(min(abs(2.0*slk), abs(sck), abs(2.0*srk)), sck)
345 else
346 as(k) = 0.
347 endif
348 endif
349 enddo
350 as(nk) = 0.
351 al(1) = s(1)
352 ar(1) = s(1)
353 al(2) = s(1)
354 do k=3,nk-1
355 al(k) = (dp(k)*s(k-1) + dp(k-1)*s(k)) * i_h12(k) &
356 + i_h0123(k)*( 2.*dp(k)*dp(k-1)*i_h12(k)*(s(k)-s(k-1)) * &
357 ( h01_h112(k) - h23_h122(k) ) &
358 + (dp(k)*as(k-1)*h23_h122(k) - dp(k-1)*as(k)*h01_h112(k)) )
359 ar(k-1) = al(k)
360 enddo
361 ar(nk-1) = s(nk)
362 al(nk) = s(nk)
363 ar(nk) = s(nk)
364 do k=2,nk-1
365 if ((pcm_layer(k)) .or. ((s(k+1)-s(k))*(s(k)-s(k-1)) <= 0.)) then
366 al(k) = s(k)
367 ar(k) = s(k)
368 else
369 da = ar(k)-al(k)
370 a6 = 6.0*s(k) - 3.0*(al(k)+ar(k))
371 if (da*a6 > da*da) then
372 al(k) = 3.0*s(k) - 2.0*ar(k)
373 elseif (da*a6 < -da*da) then
374 ar(k) = 3.0*s(k) - 2.0*al(k)
375 endif
376 endif
377 enddo
378 do k=1,nk
379 edges(k,1) = al(k)
380 edges(k,2) = ar(k)
381 enddo
382
383end subroutine hybgen_ppm_coefs
384
385!> Bound edge values by the averages of the neighboring cells.
386!!
387!! Copied from regrid_edge_values.bound_edge_values().
388subroutine bound_edge_values(N, h, u, edge_val, h_neglect, answer_date)
389 integer, intent(in) :: N !< Number of cells
390 real, dimension(N), intent(in) :: h !< Cell widths [H]
391 real, dimension(N), intent(in) :: u !< Cell averages [A]
392 real, dimension(N,2), intent(inout) :: edge_val !< Edge values [A]
393 real, intent(in) :: h_neglect !< A negligibly small width [H]
394 integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use
395
396 real :: sigma_l, sigma_c, sigma_r
397 real :: slope_x_h
398 logical :: use_2018_answers
399 integer :: k, km1, kp1
400
401 use_2018_answers = .true. ; if (present(answer_date)) use_2018_answers = (answer_date < 20190101)
402
403 do k = 1,n
404 km1 = max(1,k-1) ; kp1 = min(k+1,n)
405 slope_x_h = 0.0
406 if (use_2018_answers) then
407 sigma_l = 2.0 * ( u(k) - u(km1) ) / ( h(k) + h_neglect )
408 sigma_c = 2.0 * ( u(kp1) - u(km1) ) / ( h(km1) + 2.0*h(k) + h(kp1) + h_neglect )
409 sigma_r = 2.0 * ( u(kp1) - u(k) ) / ( h(k) + h_neglect )
410 if ( (sigma_l * sigma_r) > 0.0 ) &
411 slope_x_h = 0.5 * h(k) * sign( min(abs(sigma_l),abs(sigma_c),abs(sigma_r)), sigma_c )
412 elseif ( ((h(km1) + h(kp1)) + 2.0*h(k)) > 0.0 ) then
413 sigma_l = ( u(k) - u(km1) )
414 sigma_c = ( u(kp1) - u(km1) ) * ( h(k) / ((h(km1) + h(kp1)) + 2.0*h(k)) )
415 sigma_r = ( u(kp1) - u(k) )
416 if ( (sigma_l * sigma_r) > 0.0 ) &
417 slope_x_h = sign( min(abs(sigma_l),abs(sigma_c),abs(sigma_r)), sigma_c )
418 endif
419 if ( (u(km1)-edge_val(k,1)) * (edge_val(k,1)-u(k)) < 0.0 ) then
420 edge_val(k,1) = u(k) - sign( min( abs(slope_x_h), abs(edge_val(k,1)-u(k)) ), slope_x_h )
421 endif
422 if ( (u(kp1)-edge_val(k,2)) * (edge_val(k,2)-u(k)) < 0.0 ) then
423 edge_val(k,2) = u(k) + sign( min( abs(slope_x_h), abs(edge_val(k,2)-u(k)) ), slope_x_h )
424 endif
425 edge_val(k,1) = max( min( edge_val(k,1), max(u(km1), u(k)) ), min(u(km1), u(k)) )
426 edge_val(k,2) = max( min( edge_val(k,2), max(u(kp1), u(k)) ), min(u(kp1), u(k)) )
427 enddo
428
429end subroutine bound_edge_values
430
431!> Replace discontinuous edge values with their average when not monotonic.
432!!
433!! Copied from regrid_edge_values.check_discontinuous_edge_values().
434subroutine check_discontinuous_edge_values(N, u, edge_val)
435 integer, intent(in) :: N !< Number of cells
436 real, dimension(N), intent(in) :: u !< Cell averages [A]
437 real, dimension(N,2), intent(inout) :: edge_val !< Edge values [A]
438
439 integer :: k
440 real :: u0_avg
441
442 do k = 1,n-1
443 if ( (edge_val(k+1,1) - edge_val(k,2)) * (u(k+1) - u(k)) < 0.0 ) then
444 u0_avg = 0.5 * ( edge_val(k,2) + edge_val(k+1,1) )
445 u0_avg = max( min( u0_avg, max(u(k), u(k+1)) ), min(u(k), u(k+1)) )
446 edge_val(k,2) = u0_avg
447 edge_val(k+1,1) = u0_avg
448 endif
449 enddo
450
452
453end module recon1d_ppm_hybgen