PPM_functions.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!> Provides functions used with the Piecewise-Parabolic-Method in the vertical ALE algorithm.
6module ppm_functions
7
8! First version was created by Laurent White, June 2008.
9! Substantially re-factored January 2016.
10
11!! @todo Re-factor PPM_boundary_extrapolation to give round-off safe and
12!! optimization independent results.
13
15
16implicit none ; private
17
19
20contains
21
22!> Builds quadratic polynomials coefficients from cell mean and edge values.
23subroutine ppm_reconstruction( N, h, u, edge_values, ppoly_coef, h_neglect, answer_date)
24 integer, intent(in) :: n !< Number of cells
25 real, dimension(N), intent(in) :: h !< Cell widths [H]
26 real, dimension(N), intent(in) :: u !< Cell averages in arbitrary coordinates [A]
27 real, dimension(N,2), intent(inout) :: edge_values !< Edge values [A]
28 real, dimension(N,3), intent(inout) :: ppoly_coef !< Polynomial coefficients, mainly [A]
29 real, intent(in) :: h_neglect !< A negligibly small width [H]
30 integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use
31
32 ! Local variables
33 integer :: k ! Loop index
34 real :: edge_l, edge_r ! Edge values (left and right) [A]
35
36 ! PPM limiter
37 call ppm_limiter_standard( n, h, u, edge_values, h_neglect, answer_date=answer_date )
38
39 ! Loop over all cells
40 do k = 1,n
41
42 edge_l = edge_values(k,1)
43 edge_r = edge_values(k,2)
44
45 ! Store polynomial coefficients
46 ppoly_coef(k,1) = edge_l
47 ppoly_coef(k,2) = 4.0 * ( u(k) - edge_l ) + 2.0 * ( u(k) - edge_r )
48 ppoly_coef(k,3) = 3.0 * ( ( edge_r - u(k) ) + ( edge_l - u(k) ) )
49
50 enddo
51
52end subroutine ppm_reconstruction
53
54!> Adjusts edge values using the standard PPM limiter (Colella & Woodward, JCP 1984)
55!! after first checking that the edge values are bounded by neighbors cell averages
56!! and that the edge values are monotonic between cell averages.
57subroutine ppm_limiter_standard( N, h, u, edge_values, h_neglect, answer_date )
58 integer, intent(in) :: N !< Number of cells
59 real, dimension(:), intent(in) :: h !< cell widths (size N) [H]
60 real, dimension(:), intent(in) :: u !< cell average properties (size N) [A]
61 real, dimension(:,:), intent(inout) :: edge_values !< Potentially modified edge values [A]
62 real, intent(in) :: h_neglect !< A negligibly small width [H]
63 integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use
64
65 ! Local variables
66 integer :: k ! Loop index
67 real :: u_l, u_c, u_r ! Cell averages (left, center and right) [A]
68 real :: edge_l, edge_r ! Edge values (left and right) [A]
69 real :: expr1, expr2 ! Temporary expressions [A2]
70
71 ! Bound edge values
72 call bound_edge_values( n, h, u, edge_values, h_neglect, answer_date=answer_date )
73
74 ! Make discontinuous edge values monotonic
75 call check_discontinuous_edge_values( n, u, edge_values )
76
77 ! Loop on interior cells to apply the standard
78 ! PPM limiter (Colella & Woodward, JCP 84)
79 do k = 2,n-1
80
81 ! Get cell averages
82 u_l = u(k-1)
83 u_c = u(k)
84 u_r = u(k+1)
85
86 edge_l = edge_values(k,1)
87 edge_r = edge_values(k,2)
88
89 if ( (u_r - u_c)*(u_c - u_l) <= 0.0) then
90 ! Flatten extremum
91 edge_l = u_c
92 edge_r = u_c
93 else
94 expr1 = 3.0 * (edge_r - edge_l) * ( (u_c - edge_l) + (u_c - edge_r))
95 expr2 = (edge_r - edge_l) * (edge_r - edge_l)
96 if ( expr1 > expr2 ) then
97 ! Place extremum at right edge of cell by adjusting left edge value
98 edge_l = u_c + 2.0 * ( u_c - edge_r )
99 edge_l = max( min( edge_l, max(u_l, u_c) ), min(u_l, u_c) ) ! In case of round off
100 elseif ( expr1 < -expr2 ) then
101 ! Place extremum at left edge of cell by adjusting right edge value
102 edge_r = u_c + 2.0 * ( u_c - edge_l )
103 edge_r = max( min( edge_r, max(u_r, u_c) ), min(u_r, u_c) ) ! In case of round off
104 endif
105 endif
106 ! This checks that the difference in edge values is representable
107 ! and avoids overshoot problems due to round off.
108 !### The 1.e-60 needs to have units of [A], so this dimensionally inconsistent.
109 if ( abs( edge_r - edge_l )<max(1.e-60,epsilon(u_c)*abs(u_c)) ) then
110 edge_l = u_c
111 edge_r = u_c
112 endif
113
114 edge_values(k,1) = edge_l
115 edge_values(k,2) = edge_r
116
117 enddo ! end loop on interior cells
118
119 ! PCM within boundary cells
120 edge_values(1,:) = u(1)
121 edge_values(n,:) = u(n)
122
123end subroutine ppm_limiter_standard
124
125!> Adjusts edge values using the original monotonicity constraint (Colella & Woodward, JCP 1984)
126!! Based on hybgen_ppm_coefs
127subroutine ppm_monotonicity( N, u, edge_values )
128 integer, intent(in) :: n !< Number of cells
129 real, dimension(:), intent(in) :: u !< cell average properties (size N) [A]
130 real, dimension(:,:), intent(inout) :: edge_values !< Potentially modified edge values [A]
131
132 ! Local variables
133 integer :: k ! Loop index
134 real :: a6, da ! Normalized scalar curvature and slope [A]
135
136 ! Loop on interior cells to impose monotonicity
137 ! Eq. 1.10 of (Colella & Woodward, JCP 84)
138 do k = 2,n-1
139 if (((u(k+1)-u(k))*(u(k)-u(k-1)) <= 0.)) then !local extremum
140 edge_values(k,1) = u(k)
141 edge_values(k,2) = u(k)
142 else
143 da = edge_values(k,2)-edge_values(k,1)
144 a6 = 6.0*u(k) - 3.0*(edge_values(k,1)+edge_values(k,2))
145 if (da*a6 > da*da) then !peak in right half of zone
146 edge_values(k,1) = 3.0*u(k) - 2.0*edge_values(k,2)
147 elseif (da*a6 < -da*da) then !peak in left half of zone
148 edge_values(k,2) = 3.0*u(k) - 2.0*edge_values(k,1)
149 endif
150 endif
151 enddo ! end loop on interior cells
152
153end subroutine ppm_monotonicity
154
155!------------------------------------------------------------------------------
156!> Reconstruction by parabolas within boundary cells
157subroutine ppm_boundary_extrapolation( N, h, u, edge_values, ppoly_coef, h_neglect)
158!------------------------------------------------------------------------------
159! Reconstruction by parabolas within boundary cells.
160!
161! The following explanations apply to the left boundary cell. The same
162! reasoning holds for the right boundary cell.
163!
164! A parabola needs to be built in the cell and requires three degrees of
165! freedom, which are the right edge value and slope and the cell average.
166! The right edge values and slopes are taken to be that of the neighboring
167! cell (i.e., the left edge value and slope of the neighboring cell).
168! The resulting parabola is not necessarily monotonic and the traditional
169! PPM limiter is used to modify one of the edge values in order to yield
170! a monotonic parabola.
171!
172! N: number of cells in grid
173! h: thicknesses of grid cells
174! u: cell averages to use in constructing piecewise polynomials
175! edge_values : edge values of piecewise polynomials
176! ppoly_coef : coefficients of piecewise polynomials
177!
178! It is assumed that the size of the array 'u' is equal to the number of cells
179! defining 'grid' and 'ppoly'. No consistency check is performed here.
180!------------------------------------------------------------------------------
181
182 ! Arguments
183 integer, intent(in) :: n !< Number of cells
184 real, dimension(:), intent(in) :: h !< cell widths (size N) [H]
185 real, dimension(:), intent(in) :: u !< cell averages (size N) [A]
186 real, dimension(:,:), intent(inout) :: edge_values !< edge values of piecewise polynomials [A]
187 real, dimension(:,:), intent(inout) :: ppoly_coef !< coefficients of piecewise polynomials, mainly [A]
188 real, intent(in) :: h_neglect !< A negligibly small width for
189 !! the purpose of cell reconstructions [H]
190
191 ! Local variables
192 integer :: i0, i1
193 real :: u0, u1 ! Average concentrations in the two neighboring cells [A]
194 real :: h0, h1 ! Thicknesses of the two neighboring cells [H]
195 real :: a, b, c ! An edge value, normalized slope and normalized curvature
196 ! of a reconstructed distribution [A]
197 real :: u0_l, u0_r ! Edge values of a neighboring cell [A]
198 real :: u1_l, u1_r ! Neighboring cell slopes renormalized by the thickness of
199 ! the cell being worked on [A]
200 real :: slope ! The normalized slope [A]
201 real :: exp1, exp2 ! Temporary expressions [A2]
202
203 ! ----- Left boundary -----
204 i0 = 1
205 i1 = 2
206 h0 = h(i0)
207 h1 = h(i1)
208 u0 = u(i0)
209 u1 = u(i1)
210
211 ! Compute the left edge slope in neighboring cell and express it in
212 ! the global coordinate system
213 b = ppoly_coef(i1,2)
214 u1_r = b *((h0+h_neglect)/(h1+h_neglect)) ! derivative evaluated at xi = 0.0,
215 ! expressed w.r.t. xi (local coord. system)
216
217 ! Limit the right slope by the PLM limited slope
218 slope = 2.0 * ( u1 - u0 )
219 if ( abs(u1_r) > abs(slope) ) then
220 u1_r = slope
221 endif
222
223 ! The right edge value in the boundary cell is taken to be the left
224 ! edge value in the neighboring cell
225 u0_r = edge_values(i1,1)
226
227 ! Given the right edge value and slope, we determine the left
228 ! edge value and slope by computing the parabola as determined by
229 ! the right edge value and slope and the boundary cell average
230 u0_l = 3.0 * u0 + 0.5 * u1_r - 2.0 * u0_r
231
232 ! Apply the traditional PPM limiter
233 exp1 = (u0_r - u0_l) * (u0 - 0.5*(u0_l+u0_r))
234 exp2 = (u0_r - u0_l) * (u0_r - u0_l) / 6.0
235
236 if ( exp1 > exp2 ) then
237 u0_l = 3.0 * u0 - 2.0 * u0_r
238 endif
239
240 if ( exp1 < -exp2 ) then
241 u0_r = 3.0 * u0 - 2.0 * u0_l
242 endif
243
244 edge_values(i0,1) = u0_l
245 edge_values(i0,2) = u0_r
246
247 a = u0_l
248 b = 6.0 * u0 - 4.0 * u0_l - 2.0 * u0_r
249 c = 3.0 * ( u0_r + u0_l - 2.0 * u0 )
250
251 ppoly_coef(i0,1) = a
252 ppoly_coef(i0,2) = b
253 ppoly_coef(i0,3) = c
254
255 ! ----- Right boundary -----
256 i0 = n-1
257 i1 = n
258 h0 = h(i0)
259 h1 = h(i1)
260 u0 = u(i0)
261 u1 = u(i1)
262
263 ! Compute the right edge slope in neighboring cell and express it in
264 ! the global coordinate system
265 b = ppoly_coef(i0,2)
266 c = ppoly_coef(i0,3)
267 u1_l = (b + 2*c) ! derivative evaluated at xi = 1.0
268 u1_l = u1_l * ((h1+h_neglect)/(h0+h_neglect))
269
270 ! Limit the left slope by the PLM limited slope
271 slope = 2.0 * ( u1 - u0 )
272 if ( abs(u1_l) > abs(slope) ) then
273 u1_l = slope
274 endif
275
276 ! The left edge value in the boundary cell is taken to be the right
277 ! edge value in the neighboring cell
278 u0_l = edge_values(i0,2)
279
280 ! Given the left edge value and slope, we determine the right
281 ! edge value and slope by computing the parabola as determined by
282 ! the left edge value and slope and the boundary cell average
283 u0_r = 3.0 * u1 - 0.5 * u1_l - 2.0 * u0_l
284
285 ! Apply the traditional PPM limiter
286 exp1 = (u0_r - u0_l) * (u1 - 0.5*(u0_l+u0_r))
287 exp2 = (u0_r - u0_l) * (u0_r - u0_l) / 6.0
288
289 if ( exp1 > exp2 ) then
290 u0_l = 3.0 * u1 - 2.0 * u0_r
291 endif
292
293 if ( exp1 < -exp2 ) then
294 u0_r = 3.0 * u1 - 2.0 * u0_l
295 endif
296
297 edge_values(i1,1) = u0_l
298 edge_values(i1,2) = u0_r
299
300 a = u0_l
301 b = 6.0 * u1 - 4.0 * u0_l - 2.0 * u0_r
302 c = 3.0 * ( u0_r + u0_l - 2.0 * u1 )
303
304 ppoly_coef(i1,1) = a
305 ppoly_coef(i1,2) = b
306 ppoly_coef(i1,3) = c
307
308end subroutine ppm_boundary_extrapolation
309
310end module ppm_functions