Recon1d_MPLM_WA_poly.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!> Monotonized Piecewise Linear Method 1D reconstruction using polynomial representation
6!!
7!! This implementation of PLM follows White and Adcroft, 2008 \cite white2008.
8!! The PLM slopes are first limited following Colella and Woodward, 1984, but are then
9!! further limited to ensure the edge values moving across cell boundaries are monotone.
10!! The first and last cells are always limited to PCM.
11!!
12!! This stores and evaluates the reconstruction using a polynomial representation which is
13!! not preferred but was the form used in OM4.
15
16use recon1d_mplm_wa, only : mplm_wa, testing
17
18implicit none ; private
19
20public mplm_wa_poly, testing
21
22!> Limited Monotonic PLM reconstruction following White and Adcroft, 2008
23!!
24!! The source for the methods ultimately used by this class are:
25!! - init() *locally defined
26!! - reconstruct() *locally defined
27!! - average() *locally defined
28!! - f() -> recon1d_mplm_wa -> recon1d_plm_cw.f()
29!! - dfdx() -> recon1d_mplm_wa -> recon1d_plm_cw.dfdx()
30!! - check_reconstruction() *locally defined
31!! - unit_tests() *locally defined
32!! - destroy() -> recon1d_mplm_wa -> recon1d_plm_cw.destroy()
33!! - remap_to_sub_grid() *locally defined
34!! - init_parent() -> init()
35!! - reconstruct_parent() -> reconstruct()
36type, extends (mplm_wa) :: mplm_wa_poly
37
38 ! Legacy representation
39 integer :: degree !< Degree of polynomial used in legacy representation
40 real, allocatable, dimension(:,:) :: poly_coef !< Polynomial coefficients in legacy representation
41
42contains
43 !> Implementation of the MPLM_WA_poly initialization
44 procedure :: init => init
45 !> Implementation of the MPLM_WA_poly reconstruction
46 procedure :: reconstruct => reconstruct
47 !> Implementation of the MPLM_WA_poly average over an interval [A]
48 procedure :: average => average
49 !> Implementation of check reconstruction for the MPLM_WA_poly reconstruction
51 !> Implementation of unit tests for the MPLM_WA_poly reconstruction
52 procedure :: unit_tests => unit_tests
53
54 !> Duplicate interface to init()
55 procedure :: init_parent => init
56 !> Duplicate interface to reconstruct()
57 procedure :: reconstruct_parent => reconstruct
58
59#undef USE_BASE_CLASS_REMAP
60#ifndef USE_BASE_CLASS_REMAP
61! This block is here to test whether the compiler can do better if we have local copies of
62! the remapping functions.
63 !> Remaps the column to subgrid h_sub
65#endif
66
67end type mplm_wa_poly
68
69contains
70
71!> Initialize a 1D PLM reconstruction for n cells
72subroutine init(this, n, h_neglect, check)
73 class(mplm_wa_poly), intent(out) :: this !< This reconstruction
74 integer, intent(in) :: n !< Number of cells in this column
75 real, optional, intent(in) :: h_neglect !< A negligibly small width used in cell reconstructions [H]
76 logical, optional, intent(in) :: check !< If true, enable some consistency checking
77
78 this%n = n
79
80 allocate( this%u_mean(n) )
81 allocate( this%ul(n) )
82 allocate( this%ur(n) )
83
84 this%h_neglect = tiny( this%u_mean(1) )
85 if (present(h_neglect)) this%h_neglect = h_neglect
86 this%check = .false.
87 if (present(check)) this%check = check
88
89 this%degree = 2
90 allocate( this%poly_coef(n,2) )
91
92end subroutine init
93
94!> Calculate a 1D MPLM_WA_poly reconstructions based on h(:) and u(:)
95subroutine reconstruct(this, h, u)
96 class(mplm_wa_poly), intent(inout) :: this !< This reconstruction
97 real, intent(in) :: h(*) !< Grid spacing (thickness) [typically H]
98 real, intent(in) :: u(*) !< Cell mean values [A]
99 ! Local variables
100 real :: slp(this%n) ! The PLM slopes (difference across cell) [A]
101 real :: mslp(this%n) ! The monotonized PLM slopes [A]
102 real :: e_r, edge ! Edge values [A]
103 real :: almost_one ! A value that is slightly smaller than 1 [nondim]
104 integer :: k, n
105
106 n = this%n
107
108 ! Loop over all cells
109 do k = 1, n
110 this%u_mean(k) = u(k)
111 enddo
112
113 ! Loop on interior cells
114 do k = 2, n-1
115 slp(k) = plm_slope_wa(h(k-1), h(k), h(k+1), this%h_neglect, u(k-1), u(k), u(k+1))
116 enddo ! end loop on interior cells
117
118 ! Boundary cells use PCM. Extrapolation is handled after monotonization.
119 slp(1) = 0.
120 slp(n) = 0.
121
122 ! This loop adjusts the slope so that edge values are monotonic.
123 do k = 2, n-1
124 mslp(k) = plm_monotonized_slope( u(k-1), u(k), u(k+1), slp(k-1), slp(k), slp(k+1) )
125 enddo ! end loop on interior cells
126 mslp(1) = 0.
127 mslp(n) = 0.
128
129 ! Store and return edge values and polynomial coefficients.
130 almost_one = 1. - epsilon(e_r)
131 this%ul(1) = u(1)
132 this%ur(1) = u(1)
133 this%poly_coef(1,1) = u(1)
134 this%poly_coef(1,2) = 0.
135 do k = 2, n-1
136 this%ul(k) = u(k) - 0.5 * mslp(k) ! Left edge value of cell k
137 this%ur(k) = u(k) + 0.5 * mslp(k) ! Right edge value of cell k
138
139 this%poly_coef(k,1) = this%ul(k)
140 this%poly_coef(k,2) = this%ur(k) - this%ul(k)
141 ! Check to see if this evaluation of the polynomial at x=1 would be
142 ! monotonic w.r.t. the next cell's edge value. If not, scale back!
143 edge = this%poly_coef(k,2) + this%poly_coef(k,1)
144 e_r = u(k+1) - 0.5 * sign( mslp(k+1), slp(k+1) )
145 if ( (edge-u(k))*(e_r-edge)<0.) then
146 this%poly_coef(k,2) = this%poly_coef(k,2) * almost_one
147 endif
148 enddo
149 this%ul(n) = u(n)
150 this%ur(n) = u(n)
151 this%poly_coef(n,1) = u(n)
152 this%poly_coef(n,2) = 0.
153
154end subroutine reconstruct
155
156!> Returns a limited PLM slope following White and Adcroft, 2008, in the same arbitrary
157!! units [A] as the input values.
158!! Note that this is not the same as the Colella and Woodward method.
159real elemental pure function plm_slope_wa(h_l, h_c, h_r, h_neglect, u_l, u_c, u_r)
160 real, intent(in) :: h_l !< Thickness of left cell in arbitrary grid thickness units [H]
161 real, intent(in) :: h_c !< Thickness of center cell in arbitrary grid thickness units [H]
162 real, intent(in) :: h_r !< Thickness of right cell in arbitrary grid thickness units [H]
163 real, intent(in) :: h_neglect !< A negligible thickness [H]
164 real, intent(in) :: u_l !< Value of left cell in arbitrary units [A]
165 real, intent(in) :: u_c !< Value of center cell in arbitrary units [A]
166 real, intent(in) :: u_r !< Value of right cell in arbitrary units [A]
167 ! Local variables
168 real :: sigma_l, sigma_c, sigma_r ! Left, central and right slope estimates as
169 ! differences across the cell [A]
170 real :: u_min, u_max ! Minimum and maximum value across cell [A]
171
172 ! Side differences
173 sigma_r = u_r - u_c
174 sigma_l = u_c - u_l
175
176 ! Quasi-second order difference
177 sigma_c = 2.0 * ( u_r - u_l ) * ( h_c / ( h_l + 2.0*h_c + h_r + h_neglect) )
178
179 ! Limit slope so that reconstructions are bounded by neighbors
180 u_min = min( u_l, u_c, u_r )
181 u_max = max( u_l, u_c, u_r )
182 if ( (sigma_l * sigma_r) > 0.0 ) then
183 ! This limits the slope so that the edge values are bounded by the
184 ! two cell averages spanning the edge.
185 plm_slope_wa = sign( min( abs(sigma_c), 2.*min( u_c - u_min, u_max - u_c ) ), sigma_c )
186 else
187 ! Extrema in the mean values require a PCM reconstruction avoid generating
188 ! larger extreme values.
189 plm_slope_wa = 0.0
190 endif
191
192 ! This block tests to see if roundoff causes edge values to be out of bounds
193 if (u_c - 0.5*abs(plm_slope_wa) < u_min .or. u_c + 0.5*abs(plm_slope_wa) > u_max) then
194 plm_slope_wa = plm_slope_wa * ( 1. - epsilon(plm_slope_wa) )
195 endif
196
197 ! An attempt to avoid inconsistency when the values become unrepresentable.
198 ! ### The following 1.E-140 is dimensionally inconsistent. A newer version of
199 ! PLM is progress that will avoid the need for such rounding.
200 if (abs(plm_slope_wa) < 1.e-140) plm_slope_wa = 0.
201
202end function plm_slope_wa
203
204!> Returns a limited PLM slope following Colella and Woodward 1984, in the same
205!! arbitrary units as the input values [A].
206real elemental pure function plm_monotonized_slope(u_l, u_c, u_r, s_l, s_c, s_r)
207 real, intent(in) :: u_l !< Value of left cell in arbitrary units [A]
208 real, intent(in) :: u_c !< Value of center cell in arbitrary units [A]
209 real, intent(in) :: u_r !< Value of right cell in arbitrary units [A]
210 real, intent(in) :: s_l !< PLM slope of left cell [A]
211 real, intent(in) :: s_c !< PLM slope of center cell [A]
212 real, intent(in) :: s_r !< PLM slope of right cell [A]
213 ! Local variables
214 real :: e_r, e_l, edge ! Right, left and temporary edge values [A]
215 real :: almost_two ! The number 2, almost [nondim]
216 real :: slp ! Magnitude of PLM central slope [A]
217
218 almost_two = 2. * ( 1. - epsilon(s_c) )
219
220 ! Edge values of neighbors abutting this cell
221 e_r = u_l + 0.5*s_l
222 e_l = u_r - 0.5*s_r
223 slp = abs(s_c)
224
225 ! Check that left edge is between right edge of cell to the left and this cell mean
226 edge = u_c - 0.5 * s_c
227 if ( ( edge - e_r ) * ( u_c - edge ) < 0. ) then
228 edge = 0.5 * ( edge + e_r )
229 slp = min( slp, abs( edge - u_c ) * almost_two )
230 endif
231
232 ! Check that right edge is between left edge of cell to the right and this cell mean
233 edge = u_c + 0.5 * s_c
234 if ( ( edge - u_c ) * ( e_l - edge ) < 0. ) then
235 edge = 0.5 * ( edge + e_l )
236 slp = min( slp, abs( edge - u_c ) * almost_two )
237 endif
238
239 plm_monotonized_slope = sign( slp, s_c )
240
241end function plm_monotonized_slope
242
243!> Average between xa and xb for cell k of a 1D PLM reconstruction [A]
244!! Note: this uses the simple polynomial form a + b * x on x E (0,1)
245!! which can overshoot at x=1
246real function average(this, k, xa, xb)
247 class(mplm_wa_poly), intent(in) :: this !< This reconstruction
248 integer, intent(in) :: k !< Cell number
249 real, intent(in) :: xa !< Start of averaging interval on element (0 to 1)
250 real, intent(in) :: xb !< End of averaging interval on element (0 to 1)
251
252 average = this%poly_coef(k,1) &
253 + this%poly_coef(k,2) * 0.5 * ( xb + xa )
254
255end function average
256
257#ifndef USE_BASE_CLASS_REMAP
258! This block is needed to enable the "bounded" to test whether the compiler can do better if we have local copies of
259! the remapping functions.
260
261!> Remaps the column to subgrid h_sub
262!!
263!! It is assumed that h_sub is a perfect sub-grid of h0, meaning each h0 cell
264!! can be constructed by joining a contiguous set of h_sub cells. The integer
265!! indices isrc_start, isrc_end, isub_src provide this mapping, and are
266!! calculated in MOM_remapping
267subroutine remap_to_sub_grid(this, h0, u0, n1, h_sub, &
268 isrc_start, isrc_end, isrc_max, isub_src, &
269 u_sub, uh_sub, u02_err)
270 class(mplm_wa_poly), intent(in) :: this !< 1-D reconstruction type
271 real, intent(in) :: h0(*) !< Source grid widths (size n0) [H]
272 real, intent(in) :: u0(*) !< Source grid widths (size n0) [H]
273 integer, intent(in) :: n1 !< Number of cells in target grid
274 real, intent(in) :: h_sub(*) !< Overlapping sub-cell thicknesses, h_sub [H]
275 integer, intent(in) :: isrc_start(*) !< Index of first sub-cell within each source cell
276 integer, intent(in) :: isrc_end(*) !< Index of last sub-cell within each source cell
277 integer, intent(in) :: isrc_max(*) !< Index of thickest sub-cell within each source cell
278 integer, intent(in) :: isub_src(*) !< Index of source cell for each sub-cell
279 real, intent(out) :: u_sub(*) !< Sub-cell cell averages (size n1) [A]
280 real, intent(out) :: uh_sub(*) !< Sub-cell cell integrals (size n1) [A H]
281 real, intent(out) :: u02_err !< Integrated reconstruction error estimates [A H]
282 ! Local variables
283 integer :: i_sub ! Index of sub-cell
284 integer :: i0 ! Index into h0(1:n0), source column
285 integer :: i_max ! Used to record which sub-cell is the largest contribution of a source cell
286 real :: dh_max ! Used to record which sub-cell is the largest contribution of a source cell [H]
287 real :: xa, xb ! Non-dimensional position within a source cell (0..1) [nondim]
288 real :: dh ! The width of the sub-cell [H]
289 real :: duh ! The total amount of accumulated stuff (u*h) [A H]
290 real :: dh0_eff ! Running sum of source cell thickness [H]
291 integer :: i0_last_thick_cell, n0
292 real :: u0_min(this%n), u0_max(this%n) ! Min/max of u0 for each source cell [A]
293 real :: ul, ur ! left/right edge values of cell i0
294
295 n0 = this%n
296
297 i0_last_thick_cell = 0
298 do i0 = 1, n0
299 ul = this%ul(i0)
300 ur = this%ur(i0)
301 u0_min(i0) = min(ul, ur)
302 u0_max(i0) = max(ul, ur)
303 if (h0(i0)>0.) i0_last_thick_cell = i0
304 enddo
305
306 ! Loop over each sub-cell to calculate average/integral values within each sub-cell.
307 ! Uses: h_sub, isub_src, h0_eff
308 ! Sets: u_sub, uh_sub
309 xa = 0.
310 dh0_eff = 0.
311 u02_err = 0.
312 do i_sub = 1, n0+n1
313
314 ! Sub-cell thickness from loop above
315 dh = h_sub(i_sub)
316
317 ! Source cell
318 i0 = isub_src(i_sub)
319
320 ! Evaluate average and integral for sub-cell i_sub.
321 ! Integral is over distance dh but expressed in terms of non-dimensional
322 ! positions with source cell from xa to xb (0 <= xa <= xb <= 1).
323 dh0_eff = dh0_eff + dh ! Cumulative thickness within the source cell
324 if (h0(i0)>0.) then
325 xb = dh0_eff / h0(i0) ! This expression yields xa <= xb <= 1.0
326 xb = min(1., xb) ! This is only needed when the total target column is wider than the source column
327 u_sub(i_sub) = this%average( i0, xa, xb )
328 else ! Vanished cell
329 xb = 1.
330 u_sub(i_sub) = u0(i0)
331 endif
332 uh_sub(i_sub) = dh * u_sub(i_sub)
333 u_sub(i_sub) = max( u_sub(i_sub), u0_min(i0) )
334 u_sub(i_sub) = min( u_sub(i_sub), u0_max(i0) )
335
336 if (isub_src(i_sub+1) /= i0) then
337 ! If the next sub-cell is in a different source cell, reset the position counters
338 dh0_eff = 0.
339 xa = 0.
340 else
341 xa = xb ! Next integral will start at end of last
342 endif
343
344 enddo
345 i_sub = n0+n1+1
346 ! Sub-cell thickness from loop above
347 dh = h_sub(i_sub)
348 ! Source cell
349 i0 = isub_src(i_sub)
350
351 ! Evaluate average and integral for sub-cell i_sub.
352 ! Integral is over distance dh but expressed in terms of non-dimensional
353 ! positions with source cell from xa to xb (0 <= xa <= xb <= 1).
354 dh0_eff = dh0_eff + dh ! Cumulative thickness within the source cell
355 if (h0(i0)>0.) then
356 xb = dh0_eff / h0(i0) ! This expression yields xa <= xb <= 1.0
357 xb = min(1., xb) ! This is only needed when the total target column is wider than the source column
358 u_sub(i_sub) = this%average( i0, xa, xb )
359 else ! Vanished cell
360 xb = 1.
361 u_sub(i_sub) = u0(i0)
362 endif
363 u_sub(i_sub) = max( u_sub(i_sub), u0_min(i0) )
364 u_sub(i_sub) = min( u_sub(i_sub), u0_max(i0) )
365 uh_sub(i_sub) = dh * u_sub(i_sub)
366
367 ! Loop over each source cell substituting the integral/average for the thickest sub-cell (within
368 ! the source cell) with the residual of the source cell integral minus the other sub-cell integrals
369 ! aka a genius algorithm for accurate conservation when remapping from Robert Hallberg (\@Hallberg-NOAA).
370 ! Uses: i0_last_thick_cell, isrc_max, h_sub, isrc_start, isrc_end, uh_sub, u0, h0
371 ! Updates: uh_sub
372 do i0 = 1, i0_last_thick_cell
373 i_max = isrc_max(i0)
374 dh_max = h_sub(i_max)
375 if (dh_max > 0.) then
376 ! duh will be the sum of sub-cell integrals within the source cell except for the thickest sub-cell.
377 duh = 0.
378 do i_sub = isrc_start(i0), isrc_end(i0)
379 if (i_sub /= i_max) duh = duh + uh_sub(i_sub)
380 enddo
381 uh_sub(i_max) = u0(i0)*h0(i0) - duh
382 u02_err = u02_err + max( abs(uh_sub(i_max)), abs(u0(i0)*h0(i0)), abs(duh) )
383 endif
384 enddo
385
386 ! This should not generally be used
387 if (this%check) then
388 if ( this%check_reconstruction(h0, u0) ) stop 912 ! A debugger is required to understand why this failed
389 endif
390
391end subroutine remap_to_sub_grid
392#endif
393
394!> Checks the MPLM_WA_poly reconstruction for consistency
395logical function check_reconstruction(this, h, u)
396 class(mplm_wa_poly), intent(in) :: this !< This reconstruction
397 real, intent(in) :: h(*) !< Grid spacing (thickness) [typically H]
398 real, intent(in) :: u(*) !< Cell mean values [A]
399 ! Local variables
400 integer :: k
401
402 check_reconstruction = .false.
403
404 do k = 1, this%n
405 if ( abs( this%u_mean(k) - u(k) ) > 0. ) check_reconstruction = .true.
406 enddo
407
408 ! Check the cell reconstruction is monotonic within each cell (it should be as a straight line)
409 do k = 1, this%n
410 if ( ( this%u_mean(k) - this%ul(k) ) * ( this%ur(k) - this%u_mean(k) ) < 0. ) check_reconstruction = .true.
411 enddo
412
413 ! Check the cell is a straight line (to within machine precision)
414 do k = 1, this%n
415 if ( abs(2. * this%u_mean(k) - ( this%ul(k) + this%ur(k) )) > epsilon(this%u_mean(1)) * &
416 max(abs(2. * this%u_mean(k)), abs(this%ul(k)), abs(this%ur(k))) ) check_reconstruction = .true.
417 enddo
418
419 ! Check bounding of right edges, w.r.t. the cell means
420 do k = 1, this%n-1
421 if ( ( this%ur(k) - this%u_mean(k) ) * ( this%u_mean(k+1) - this%ur(k) ) < 0. ) check_reconstruction = .true.
422 enddo
423
424 ! Check bounding of left edges, w.r.t. the cell means
425 do k = 2, this%n
426 if ( ( this%u_mean(k) - this%ul(k) ) * ( this%ul(k) - this%u_mean(k-1) ) < 0. ) check_reconstruction = .true.
427 enddo
428
429 ! Check bounding of left edges, w.r.t. this cell mean and the previous cell right edge
430 ! Check order of u, ur, ul
431 ! Note that in OM4 implementation, we were not consistent for top and bottom layers due
432 ! extrapolation using cell means rather than edge values
433 do k = 2, this%n-2
434 if ( ( this%ur(k) - this%u_mean(k) ) * ( this%ul(k+1) - this%ur(k) ) < 0. ) check_reconstruction = .true.
435 enddo
436
437end function check_reconstruction
438
439!> Runs PLM reconstruction unit tests and returns True for any fails, False otherwise
440logical function unit_tests(this, verbose, stdout, stderr)
441 class(mplm_wa_poly), intent(inout) :: this !< This reconstruction
442 logical, intent(in) :: verbose !< True, if verbose
443 integer, intent(in) :: stdout !< I/O channel for stdout
444 integer, intent(in) :: stderr !< I/O channel for stderr
445 ! Local variables
446 real, allocatable :: ul(:), ur(:), um(:) ! test values [A]
447 real, allocatable :: ull(:), urr(:) ! test values [A]
448 type(testing) :: test ! convenience functions
449 integer :: k
450
451 call test%set( stdout=stdout ) ! Sets the stdout channel in test
452 call test%set( stderr=stderr ) ! Sets the stderr channel in test
453 call test%set( verbose=verbose ) ! Sets the verbosity flag in test
454
455 call this%init(3)
456 call test%test( this%n /= 3, 'Setting number of levels')
457 allocate( um(3), ul(3), ur(3), ull(3), urr(3) )
458
459 call this%reconstruct( (/2.,2.,2./), (/1.,3.,5./) )
460 call test%real_arr(3, this%u_mean, (/1.,3.,5./), 'Setting cell values')
461
462 do k = 1, 3
463 ul(k) = this%f(k, 0.)
464 um(k) = this%f(k, 0.5)
465 ur(k) = this%f(k, 1.)
466 enddo
467 call test%real_arr(3, ul, (/1.,2.,5./), 'Evaluation on left edge')
468 call test%real_arr(3, um, (/1.,3.,5./), 'Evaluation in center')
469 call test%real_arr(3, ur, (/1.,4.,5./), 'Evaluation on right edge')
470
471 do k = 1, 3
472 ul(k) = this%dfdx(k, 0.)
473 um(k) = this%dfdx(k, 0.5)
474 ur(k) = this%dfdx(k, 1.)
475 enddo
476 call test%real_arr(3, ul, (/0.,2.,0./), 'dfdx on left edge')
477 call test%real_arr(3, um, (/0.,2.,0./), 'dfdx in center')
478 call test%real_arr(3, ur, (/0.,2.,0./), 'dfdx on right edge')
479
480 do k = 1, 3
481 um(k) = this%average(k, 0.5, 0.75) ! Average from x=0.25 to 0.75 in each cell
482 enddo
483 call test%real_arr(3, um, (/1.,3.25,5./), 'Return interval average')
484
485 unit_tests = test%summarize('MPLM_WA_poly:unit_tests')
486
487end function unit_tests
488
489!> \namespace recon1d_mplm_wa_poly
490!!
491
492end module recon1d_mplm_wa_poly