Recon1d_MPLM_WA.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
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 differs from recon1d_mplm_wa_poly in the internally not polynomial representations
13!! are referred to.
14module recon1d_mplm_wa
15
16use recon1d_plm_cw, only : plm_cw, testing
17
18implicit none ; private
19
20public mplm_wa, 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() -> recon1d_plm_cw.init()
26!! - reconstruct() *locally defined
27!! - average() -> recon1d_plm_cw.average()
28!! - f() -> recon1d_plm_cw.f()
29!! - dfdx() -> recon1d_plm_cw.dfdx()
30!! - check_reconstruction() *locally defined
31!! - unit_tests() *locally defined
32!! - destroy() -> recon1d_plm_cw.destroy()
33!! - remap_to_sub_grid() -> recon1d_plm_cw -> recon1d_type.remap_to_sub_grid()
34!! - init_parent() -> recon1d_plm_cw.init()
35!! - reconstruct_parent() -> reconstruct()
37
38contains
39 !> Implementation of the MPLM_WA reconstruction
40 procedure :: reconstruct => reconstruct
41 !> Implementation of check reconstruction for the MPLM_WA reconstruction
42 procedure :: check_reconstruction => check_reconstruction
43 !> Implementation of unit tests for the MPLM_WA reconstruction
44 procedure :: unit_tests => unit_tests
45
46 !> Duplicate interface to reconstruct()
47 procedure :: reconstruct_parent => reconstruct
48
49end type mplm_wa
50
51contains
52
53!> Calculate a 1D PLM reconstructions based on h(:) and u(:)
54subroutine reconstruct(this, h, u)
55 class(mplm_wa), intent(inout) :: this !< This reconstruction
56 real, intent(in) :: h(*) !< Grid spacing (thickness) [typically H]
57 real, intent(in) :: u(*) !< Cell mean values [A]
58 ! Local variables
59 real :: slp(this%n) ! The PLM slopes (difference across cell) [A]
60 real :: mslp(this%n) ! The monotonized PLM slopes [A]
61 integer :: k, n
62 real :: u_tmp, u_min, u_max ! Working values of cells [A]
63
64 n = this%n
65
66 ! Loop over all cells
67 do k = 1, n
68 this%u_mean(k) = u(k)
69 enddo
70
71 ! Loop on interior cells
72 do k = 2, n-1
73 slp(k) = plm_slope_wa(h(k-1), h(k), h(k+1), this%h_neglect, u(k-1), u(k), u(k+1))
74 enddo ! end loop on interior cells
75
76 ! Boundary cells use PCM. Extrapolation is handled after monotonization.
77 slp(1) = 0.
78 slp(n) = 0.
79
80 ! This loop adjusts the slope so that edge values are monotonic.
81 do k = 2, n-1
82 mslp(k) = plm_monotonized_slope( u(k-1), u(k), u(k+1), slp(k-1), slp(k), slp(k+1) )
83 enddo ! end loop on interior cells
84 mslp(1) = 0.
85 mslp(n) = 0.
86
87 ! Store edge values
88 this%ul(1) = u(1)
89 this%ur(1) = u(1)
90 do k = 2, n-1
91 u_tmp = u(k-1) + 0.5 * mslp(k-1) ! Right edge value of cell k-1
92 u_min = min( u(k), u_tmp )
93 u_max = max( u(k), u_tmp )
94 u_tmp = u(k) - 0.5 * mslp(k) ! Left edge value of cell k
95 this%ul(k) = max( min( u_tmp, u_max), u_min ) ! Bounded to handle roundoff
96 u_tmp = u(k+1) - 0.5 * mslp(k-1) ! Left edge value of cell k+1
97 u_min = min( u(k), u_tmp )
98 u_max = max( u(k), u_tmp )
99 u_tmp = u(k) + 0.5 * mslp(k) ! Right edge value of cell k
100 this%ur(k) = max( min( u_tmp, u_max), u_min ) ! Bounded to handle roundoff
101 enddo
102 this%ul(n) = u(n)
103 this%ur(n) = u(n)
104
105end subroutine reconstruct
106
107!> Returns a limited PLM slope following White and Adcroft, 2008, in the same arbitrary
108!! units [A] as the input values.
109!! Note that this is not the same as the Colella and Woodward method.
110real elemental pure function plm_slope_wa(h_l, h_c, h_r, h_neglect, u_l, u_c, u_r)
111 real, intent(in) :: h_l !< Thickness of left cell in arbitrary grid thickness units [H]
112 real, intent(in) :: h_c !< Thickness of center cell in arbitrary grid thickness units [H]
113 real, intent(in) :: h_r !< Thickness of right cell in arbitrary grid thickness units [H]
114 real, intent(in) :: h_neglect !< A negligible thickness [H]
115 real, intent(in) :: u_l !< Value of left cell in arbitrary units [A]
116 real, intent(in) :: u_c !< Value of center cell in arbitrary units [A]
117 real, intent(in) :: u_r !< Value of right cell in arbitrary units [A]
118 ! Local variables
119 real :: sigma_l, sigma_c, sigma_r ! Left, central and right slope estimates as
120 ! differences across the cell [A]
121 real :: u_min, u_max ! Minimum and maximum value across cell [A]
122
123 ! Side differences
124 sigma_r = u_r - u_c
125 sigma_l = u_c - u_l
126
127 ! Quasi-second order difference
128 sigma_c = 2.0 * ( u_r - u_l ) * ( h_c / ( h_l + 2.0*h_c + h_r + h_neglect) )
129
130 ! Limit slope so that reconstructions are bounded by neighbors
131 u_min = min( u_l, u_c, u_r )
132 u_max = max( u_l, u_c, u_r )
133 if ( (sigma_l * sigma_r) > 0.0 ) then
134 ! This limits the slope so that the edge values are bounded by the
135 ! two cell averages spanning the edge.
136 plm_slope_wa = sign( min( abs(sigma_c), 2.*min( u_c - u_min, u_max - u_c ) ), sigma_c )
137 else
138 ! Extrema in the mean values require a PCM reconstruction avoid generating
139 ! larger extreme values.
140 plm_slope_wa = 0.0
141 endif
142
143end function plm_slope_wa
144
145!> Returns a limited PLM slope following Colella and Woodward 1984, in the same
146!! arbitrary units as the input values [A].
147real elemental pure function plm_monotonized_slope(u_l, u_c, u_r, s_l, s_c, s_r)
148 real, intent(in) :: u_l !< Value of left cell in arbitrary units [A]
149 real, intent(in) :: u_c !< Value of center cell in arbitrary units [A]
150 real, intent(in) :: u_r !< Value of right cell in arbitrary units [A]
151 real, intent(in) :: s_l !< PLM slope of left cell [A]
152 real, intent(in) :: s_c !< PLM slope of center cell [A]
153 real, intent(in) :: s_r !< PLM slope of right cell [A]
154 ! Local variables
155 real :: neighbor_edge ! Edge value of nieghbor cell [A]
156 real :: this_edge ! Edge value of this cell [A]
157 real :: slp ! Magnitude of PLM central slope [A]
158
159 ! Comparison are made assuming +ve slopes
160 slp = abs(s_c)
161
162 ! Check that left edge is between right edge of cell to the left and this cell mean
163 neighbor_edge = u_l + 0.5 * s_l
164 this_edge = u_c - 0.5 * s_c
165 if ( ( this_edge - neighbor_edge ) * ( u_c - this_edge ) < 0. ) then
166 ! Using the midpoint works because the neighbor is similarly adjusted
167 this_edge = 0.5 * ( this_edge + neighbor_edge )
168 slp = min( slp, abs( this_edge - u_c ) * 2. )
169 endif
170
171 ! Check that right edge is between left edge of cell to the right and this cell mean
172 neighbor_edge = u_r - 0.5 * s_r
173 this_edge = u_c + 0.5 * s_c
174 if ( ( this_edge - u_c ) * ( neighbor_edge - this_edge ) < 0. ) then
175 ! Using the midpoint works because the neighbor is similarly adjusted
176 this_edge = 0.5 * ( this_edge + neighbor_edge )
177 slp = min( slp, abs( this_edge - u_c ) * 2. )
178 endif
179
180 plm_monotonized_slope = sign( slp, s_c )
181
182end function plm_monotonized_slope
183
184!> Checks the MPLM_WA reconstruction for consistency
185logical function check_reconstruction(this, h, u)
186 class(mplm_wa), intent(in) :: this !< This reconstruction
187 real, intent(in) :: h(*) !< Grid spacing (thickness) [typically H]
188 real, intent(in) :: u(*) !< Cell mean values [A]
189 ! Local variables
190 integer :: k
191
192 check_reconstruction = .false.
193
194 do k = 1, this%n
195 if ( abs( this%u_mean(k) - u(k) ) > 0. ) check_reconstruction = .true.
196 enddo
197
198 ! Check the cell reconstruction is monotonic within each cell (it should be as a straight line)
199 do k = 1, this%n
200 if ( ( this%u_mean(k) - this%ul(k) ) * ( this%ur(k) - this%u_mean(k) ) < 0. ) check_reconstruction = .true.
201 enddo
202
203 ! This next test fails abysmally!
204 ! Using intel/2023.2.0 on gaea, MOM_remapping:test_recon_consistency iter=6
205 ! um~0.581492556923472 ul~0.402083491713151 ur~0.749082615698503
206 ! Check the cell is a straight line (to within machine precision)
207! do k = 1, this%n
208! if ( abs(2. * this%u_mean(k) - ( this%ul(k) + this%ur(k) )) > epsilon(this%u_mean(1)) * &
209! max(abs(2. * this%u_mean(k)), abs(this%ul(k)), abs(this%ur(k))) ) check_reconstruction = .true.
210! enddo
211
212 ! Check bounding of right edges, w.r.t. the cell means
213 do k = 1, this%n-1
214 if ( ( this%ur(k) - this%u_mean(k) ) * ( this%u_mean(k+1) - this%ur(k) ) < 0. ) check_reconstruction = .true.
215 enddo
216
217 ! Check bounding of left edges, w.r.t. the cell means
218 do k = 2, this%n
219 if ( ( this%u_mean(k) - this%ul(k) ) * ( this%ul(k) - this%u_mean(k-1) ) < 0. ) check_reconstruction = .true.
220 enddo
221
222 ! Check order of u, ur, ul
223 do k = 1, this%n-1
224 if ( ( this%ur(k) - this%u_mean(k) ) * ( this%ul(k+1) - this%ur(k) ) < 0. ) check_reconstruction = .true.
225 enddo
226
227 ! Check bounding of left edges, w.r.t. this cell mean and the previous cell right edge
228 do k = 2, this%n
229 if ( ( this%u_mean(k) - this%ul(k) ) * ( this%ul(k) - this%ur(k-1) ) < 0. ) check_reconstruction = .true.
230 enddo
231
232end function check_reconstruction
233
234!> Runs PLM reconstruction unit tests and returns True for any fails, False otherwise
235logical function unit_tests(this, verbose, stdout, stderr)
236 class(mplm_wa), intent(inout) :: this !< This reconstruction
237 logical, intent(in) :: verbose !< True, if verbose
238 integer, intent(in) :: stdout !< I/O channel for stdout
239 integer, intent(in) :: stderr !< I/O channel for stderr
240 ! Local variables
241 real, allocatable :: ul(:), ur(:), um(:) ! test values [A]
242 real, allocatable :: ull(:), urr(:) ! test values [A]
243 type(testing) :: test ! convenience functions
244 integer :: k
245
246 call test%set( stdout=stdout ) ! Sets the stdout channel in test
247 call test%set( stderr=stderr ) ! Sets the stderr channel in test
248 call test%set( verbose=verbose ) ! Sets the verbosity flag in test
249
250 call this%init(3)
251 call test%test( this%n /= 3, 'Setting number of levels')
252 allocate( um(3), ul(3), ur(3), ull(3), urr(3) )
253
254 call this%reconstruct( (/2.,2.,2./), (/1.,3.,5./) )
255 call test%real_arr(3, this%u_mean, (/1.,3.,5./), 'Setting cell values')
256
257 do k = 1, 3
258 ul(k) = this%f(k, 0.)
259 um(k) = this%f(k, 0.5)
260 ur(k) = this%f(k, 1.)
261 enddo
262 call test%real_arr(3, ul, (/1.,2.,5./), 'Evaluation on left edge')
263 call test%real_arr(3, um, (/1.,3.,5./), 'Evaluation in center')
264 call test%real_arr(3, ur, (/1.,4.,5./), 'Evaluation on right edge')
265
266 do k = 1, 3
267 ul(k) = this%dfdx(k, 0.)
268 um(k) = this%dfdx(k, 0.5)
269 ur(k) = this%dfdx(k, 1.)
270 enddo
271 call test%real_arr(3, ul, (/0.,2.,0./), 'dfdx on left edge')
272 call test%real_arr(3, um, (/0.,2.,0./), 'dfdx in center')
273 call test%real_arr(3, ur, (/0.,2.,0./), 'dfdx on right edge')
274
275 do k = 1, 3
276 um(k) = this%average(k, 0.5, 0.75) ! Average from x=0.25 to 0.75 in each cell
277 enddo
278 call test%real_arr(3, um, (/1.,3.25,5./), 'Return interval average')
279
280 unit_tests = test%summarize('MPLM_WA:unit_tests')
281
282end function unit_tests
283
284!> \namespace recon1d_mplm_wa
285!!
286
287end module recon1d_mplm_wa