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