Recon1d_MPLM_CWK.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 Linear Method 1D reconstruction in index space
6!!
7!! This implementation of PLM follows Colella and Woodward, 1984 \cite colella1984, except for assuming
8!! uniform resolution so that the method is independent of grid spacing. The cell-wise reconstructions
9!! are limited so that the edge values (which are also the extrema in a cell) are bounded by the neighbors.
10!! The first and last cells are always limited to PCM.
12
13use recon1d_type, only : testing
14use recon1d_plm_cwk, only : plm_cwk
15
16implicit none ; private
17
18public mplm_cwk, testing
19
20!> PLM reconstruction following Colella and Woodward, 1984
21!!
22!! Implemented by extending recon1d_plm_cwk.
23!!
24!! The source for the methods ultimately used by this class are:
25!! - init() -> recon1d_plm_cwk -> recon1d_plm_cw.init()
26!! - reconstruct() *locally defined
27!! - average() -> recon1d_plm_cwk -> recon1d_plm_cw.average()
28!! - f() -> recon1d_plm_cwk -> recon1d_plm_cw.f()
29!! - dfdx() -> recon1d_plm_cwk -> recon1d_plm_cw.dfdx()
30!! - check_reconstruction() *locally defined
31!! - unit_tests() *locally defined
32!! - destroy() -> recon1d_plm_cwk -> recon1d_plm_cw.destroy()
33!! - remap_to_sub_grid() -> recon1d_type.remap_to_sub_grid()
34!! - init_parent() -> init()
35!! - reconstruct_parent() -> reconstruct()
36type, extends (plm_cwk) :: mplm_cwk
37
38contains
39 !> Implementation of the MPLM_CWK reconstruction
40 procedure :: reconstruct => reconstruct
41 !> Implementation of check reconstruction for the MPLM_CWK reconstruction
43 !> Implementation of unit tests for the MPLM_CWK reconstruction
44 procedure :: unit_tests => unit_tests
45
46 !> Duplicate interface to reconstruct()
47 procedure :: reconstruct_parent => reconstruct
48end type mplm_cwk
49
50contains
51
52!> Calculate a 1D PLM reconstructions based on h(:) and u(:)
53subroutine reconstruct(this, h, u)
54 class(mplm_cwk), intent(inout) :: this !< This reconstruction
55 real, intent(in) :: h(*) !< Grid spacing (thickness) [typically H]
56 real, intent(in) :: u(*) !< Cell mean values [A]
57 ! Local variables
58 real :: slp ! The PLM slopes (difference across cell) [A]
59 real :: sigma_l, sigma_c, sigma_r ! Left, central and right slope estimates as
60 ! differences across the cell [A]
61 real :: u_min, u_max ! Minimum and maximum value across cell [A]
62 real :: u_l, u_r, u_c ! Left, right, and center values [A]
63 real :: u_e(this%n+1) ! Average of edge values [A]
64 integer :: k, n
65
66 n = this%n
67
68 ! Loop over all cells
69 do k = 1, n
70 this%u_mean(k) = u(k)
71 enddo
72
73 ! Boundary cells use PCM
74 this%ul(1) = u(1)
75 this%ur(1) = u(1)
76
77 ! Loop over interior cells
78 do k = 2, n-1
79 u_l = u(k-1)
80 u_c = u(k)
81 u_r = u(k+1)
82
83 ! Side differences
84 sigma_r = u_r - u_c
85 sigma_l = u_c - u_l
86
87 ! This is the second order slope given by equation 1.7 of
88 ! Piecewise Parabolic Method, Colella and Woodward (1984),
89 ! http://dx.doi.org/10.1016/0021-991(84)90143-8.
90 ! For uniform resolution it simplifies to ( u_r - u_l )/2 .
91 sigma_c = 0.5 * ( u_r - u_l )
92
93 ! Limit slope so that reconstructions are bounded by neighbors
94 u_min = min( u_l, u_c, u_r )
95 u_max = max( u_l, u_c, u_r )
96
97 if ( (sigma_l * sigma_r) > 0.0 ) then
98 ! This limits the slope so that the edge values are bounded by the two cell averages spanning the edge
99 slp = sign( min( abs(sigma_c), 2.*min( u_c - u_min, u_max - u_c ) ), sigma_c )
100 else
101 ! Extrema in the mean values require a PCM reconstruction
102 slp = 0.0
103 endif
104
105 ! Left edge
106 u_min = min( u_c, u_l )
107 u_max = max( u_c, u_l )
108 u_l = u_c - 0.5 * slp
109 this%ul(k) = max( min( u_l, u_max), u_min )
110
111 ! Right edge
112 u_min = min( u_c, u_r )
113 u_max = max( u_c, u_r )
114 u_r = u_c + 0.5 * slp
115 this%ur(k) = max( min( u_r, u_max), u_min )
116 enddo
117
118 ! Boundary cells use PCM
119 this%ul(n) = u(n)
120 this%ur(n) = u(n)
121
122 ! Average edge values
123 u_e(1) = this%ul(1)
124 do k = 2, n
125 u_e(k) = 0.5 * ( this%ur(k-1) + this%ul(k) )
126 enddo
127 u_e(n+1) = this%ur(n)
128
129 ! Loop over interior cells, redo PLM slope limiting using average edge as neighbor cell values
130 do k = 2, n-1
131 u_l = u_e(k)
132 u_c = u(k)
133 u_r = u_e(k+1)
134
135 ! Side differences
136 sigma_r = u_r - u_c
137 sigma_l = u_c - u_l
138
139 ! This is the second order slope given by equation 1.7 of
140 ! Piecewise Parabolic Method, Colella and Woodward (1984),
141 ! http://dx.doi.org/10.1016/0021-991(84)90143-8.
142 ! For uniform resolution it simplifies to ( u_r - u_l )/2 .
143 sigma_c = this%ur(k) - this%ul(k)
144
145 ! Limit slope so that reconstructions are bounded by neighbors
146 u_min = min( u_l, u_c, u_r )
147 u_max = max( u_l, u_c, u_r )
148
149 if ( (sigma_l * sigma_r) > 0.0 ) then
150 ! This limits the slope so that the edge values are bounded by the two cell averages spanning the edge
151 slp = sign( min( abs(sigma_c), 2.*min( u_c - u_min, u_max - u_c ) ), sigma_c )
152 else
153 ! Extrema in the mean values require a PCM reconstruction
154 slp = 0.0
155 endif
156
157 ! Left edge
158 u_min = min( u_c, u_l )
159 u_max = max( u_c, u_l )
160 u_l = u_c - 0.5 * slp
161 this%ul(k) = max( min( u_l, u_max), u_min )
162
163 ! Right edge
164 u_min = min( u_c, u_r )
165 u_max = max( u_c, u_r )
166 u_r = u_c + 0.5 * slp
167 this%ur(k) = max( min( u_r, u_max), u_min )
168 enddo
169
170end subroutine reconstruct
171
172!> Checks the MPLM_CWK reconstruction for consistency
173logical function check_reconstruction(this, h, u)
174 class(mplm_cwk), intent(in) :: this !< This reconstruction
175 real, intent(in) :: h(*) !< Grid spacing (thickness) [typically H]
176 real, intent(in) :: u(*) !< Cell mean values [A]
177 ! Local variables
178 integer :: k
179
180 check_reconstruction = .false.
181
182 do k = 1, this%n
183 if ( abs( this%u_mean(k) - u(k) ) > 0. ) check_reconstruction = .true.
184 enddo
185
186 ! Check the cell reconstruction is monotonic within each cell (it should be as a straight line)
187 do k = 1, this%n
188 if ( ( this%u_mean(k) - this%ul(k) ) * ( this%ur(k) - this%u_mean(k) ) < 0. ) check_reconstruction = .true.
189 enddo
190
191 ! Check the cell is a straight line (to within machine precision)
192 do k = 1, this%n
193 if ( abs(2. * this%u_mean(k) - ( this%ul(k) + this%ur(k) )) > epsilon(this%u_mean(1)) * &
194 max(abs(2. * this%u_mean(k)), abs(this%ul(k)), abs(this%ur(k))) ) check_reconstruction = .true.
195 enddo
196
197 ! Check bounding of right edges, w.r.t. the cell means
198 do k = 1, this%n-1
199 if ( ( this%ur(k) - this%u_mean(k) ) * ( this%u_mean(k+1) - this%ur(k) ) < 0. ) check_reconstruction = .true.
200 enddo
201
202 ! Check bounding of left edges, w.r.t. the cell means
203 do k = 2, this%n
204 if ( ( this%u_mean(k) - this%ul(k) ) * ( this%ul(k) - this%u_mean(k-1) ) < 0. ) check_reconstruction = .true.
205 enddo
206
207 ! Check bounding of right edges, w.r.t. this cell mean and the next cell left edge
208 do k = 1, this%n-1
209 if ( ( this%ur(k) - this%u_mean(k) ) * ( this%ul(k+1) - this%ur(k) ) < 0. ) check_reconstruction = .true.
210 enddo
211
212 ! Check bounding of left edges, w.r.t. this cell mean and the previous cell right edge
213 do k = 2, this%n
214 if ( ( this%u_mean(k) - this%ul(k) ) * ( this%ul(k) - this%ur(k-1) ) < 0. ) check_reconstruction = .true.
215 enddo
216
217end function check_reconstruction
218
219!> Runs PLM reconstruction unit tests and returns True for any fails, False otherwise
220logical function unit_tests(this, verbose, stdout, stderr)
221 class(mplm_cwk), intent(inout) :: this !< This reconstruction
222 logical, intent(in) :: verbose !< True, if verbose
223 integer, intent(in) :: stdout !< I/O channel for stdout
224 integer, intent(in) :: stderr !< I/O channel for stderr
225 ! Local variables
226 real, allocatable :: ul(:), ur(:), um(:) ! test values [A]
227 real, allocatable :: ull(:), urr(:) ! test values [A]
228 type(testing) :: test ! convenience functions
229 integer :: k
230
231 call test%set( stdout=stdout ) ! Sets the stdout channel in test
232 call test%set( stderr=stderr ) ! Sets the stderr channel in test
233 call test%set( verbose=verbose ) ! Sets the verbosity flag in test
234
235 call this%init(3)
236 call test%test( this%n /= 3, 'Setting number of levels')
237 allocate( um(3), ul(3), ur(3), ull(3), urr(3) )
238
239 call this%reconstruct( (/2.,2.,2./), (/1.,3.,5./) )
240 call test%real_arr(3, this%u_mean, (/1.,3.,5./), 'Setting cell values')
241
242 do k = 1, 3
243 ul(k) = this%f(k, 0.)
244 um(k) = this%f(k, 0.5)
245 ur(k) = this%f(k, 1.)
246 enddo
247 call test%real_arr(3, ul, (/1.,2.,5./), 'Evaluation on left edge')
248 call test%real_arr(3, um, (/1.,3.,5./), 'Evaluation in center')
249 call test%real_arr(3, ur, (/1.,4.,5./), 'Evaluation on right edge')
250
251 do k = 1, 3
252 ul(k) = this%dfdx(k, 0.)
253 um(k) = this%dfdx(k, 0.5)
254 ur(k) = this%dfdx(k, 1.)
255 enddo
256 call test%real_arr(3, ul, (/0.,2.,0./), 'dfdx on left edge')
257 call test%real_arr(3, um, (/0.,2.,0./), 'dfdx in center')
258 call test%real_arr(3, ur, (/0.,2.,0./), 'dfdx on right edge')
259
260 do k = 1, 3
261 um(k) = this%average(k, 0.5, 0.75) ! Average from x=0.25 to 0.75 in each cell
262 enddo
263 call test%real_arr(3, um, (/1.,3.25,5./), 'Return interval average')
264
265 call this%destroy()
266 deallocate( um, ul, ur, ull, urr )
267
268 allocate( um(4), ul(4), ur(4) )
269 call this%init(4)
270
271 ! These values lead to non-monotonic reconstuctions which are
272 ! valid for transport problems but not always appropriate for
273 ! remapping to arbitrary resolution grids.
274 ! The O(h^2) slopes are -, 2, 2, - and the limited
275 ! slopes are 0, 1, 1, 0 so the everywhere the reconstructions
276 ! are bounded by neighbors but ur(2) and ul(3) are out-of-order.
277 call this%reconstruct( (/1.,1.,1.,1./), (/0.,3.,4.,7./) )
278 do k = 1, 4
279 ul(k) = this%f(k, 0.)
280 ur(k) = this%f(k, 1.)
281 enddo
282 call test%real_arr(4, ul, (/0.,2.5,3.5,7./), 'Evaluation on left edge')
283 call test%real_arr(4, ur, (/0.,3.5,4.5,7./), 'Evaluation on right edge')
284
285 deallocate( um, ul, ur )
286
287 unit_tests = test%summarize('MPLM_CWK:unit_tests')
288
289end function unit_tests
290
291!> \namespace recon1d_mplm_cwk
292!!
293
294end module recon1d_mplm_cwk