MOM6
ppm_functions Module Reference

Detailed Description

Provides functions used with the Piecewise-Parabolic-Method in the vertical ALE algorithm.

Functions/Subroutines

subroutine, public ppm_reconstruction (N, h, u, ppoly_E, ppoly_coefficients)
 Builds quadratic polynomials coefficients from cell mean and edge values. More...
 
subroutine ppm_limiter_standard (N, h, u, ppoly_E)
 Adjusts edge values using the standard PPM limiter (Colella & Woodward, JCP 1984) after first checking that the edge values are bounded by neighbors cell averages and that the edge values are monotonic between cell averages. More...
 
subroutine, public ppm_boundary_extrapolation (N, h, u, ppoly_E, ppoly_coefficients)
 

Variables

real, parameter h_neglect = 1.E-30
 A tiny width that is so small that adding it to cell widths does not change the value due to a computational representation. It is used to avoid division by zero. More...
 

Function/Subroutine Documentation

◆ ppm_boundary_extrapolation()

subroutine, public ppm_functions::ppm_boundary_extrapolation ( integer, intent(in)  N,
real, dimension(:), intent(in)  h,
real, dimension(:), intent(in)  u,
real, dimension(:,:), intent(inout)  ppoly_E,
real, dimension(:,:), intent(inout)  ppoly_coefficients 
)

Definition at line 128 of file PPM_functions.F90.

References h_neglect.

Referenced by mom_remapping::build_reconstructions_1d(), mom_ale::pressure_gradient_ppm(), regrid_interp::regridding_set_ppolys(), and mom_remapping::remapping_unit_tests().

128 !------------------------------------------------------------------------------
129 ! Reconstruction by parabolas within boundary cells.
130 !
131 ! The following explanations apply to the left boundary cell. The same
132 ! reasoning holds for the right boundary cell.
133 !
134 ! A parabola needs to be built in the cell and requires three degrees of
135 ! freedom, which are the right edge value and slope and the cell average.
136 ! The right edge values and slopes are taken to be that of the neighboring
137 ! cell (i.e., the left edge value and slope of the neighboring cell).
138 ! The resulting parabola is not necessarily monotonic and the traditional
139 ! PPM limiter is used to modify one of the edge values in order to yield
140 ! a monotonic parabola.
141 !
142 ! N: number of cells in grid
143 ! h: thicknesses of grid cells
144 ! u: cell averages to use in constructing piecewise polynomials
145 ! ppoly_E : edge values of piecewise polynomials
146 ! ppoly_coefficients : coefficients of piecewise polynomials
147 !
148 ! It is assumed that the size of the array 'u' is equal to the number of cells
149 ! defining 'grid' and 'ppoly'. No consistency check is performed here.
150 !------------------------------------------------------------------------------
151 
152  ! Arguments
153  integer, intent(in) :: n ! Number of cells
154  real, dimension(:), intent(in) :: h ! cell widths (size N)
155  real, dimension(:), intent(in) :: u ! cell averages (size N)
156  real, dimension(:,:), intent(inout) :: ppoly_e
157  real, dimension(:,:), intent(inout) :: ppoly_coefficients
158 
159  ! Local variables
160  integer :: i0, i1
161  real :: u0, u1
162  real :: h0, h1
163  real :: a, b, c
164  real :: u0_l, u0_r
165  real :: u1_l, u1_r
166  real :: slope
167  real :: exp1, exp2
168 
169  ! ----- Left boundary -----
170  i0 = 1
171  i1 = 2
172  h0 = h(i0)
173  h1 = h(i1)
174  u0 = u(i0)
175  u1 = u(i1)
176 
177  ! Compute the left edge slope in neighboring cell and express it in
178  ! the global coordinate system
179  b = ppoly_coefficients(i1,2)
180  u1_r = b *((h0+h_neglect)/(h1+h_neglect)) ! derivative evaluated at xi = 0.0,
181  ! expressed w.r.t. xi (local coord. system)
182 
183  ! Limit the right slope by the PLM limited slope
184  slope = 2.0 * ( u1 - u0 )
185  if ( abs(u1_r) .GT. abs(slope) ) then
186  u1_r = slope
187  end if
188 
189  ! The right edge value in the boundary cell is taken to be the left
190  ! edge value in the neighboring cell
191  u0_r = ppoly_e(i1,1)
192 
193  ! Given the right edge value and slope, we determine the left
194  ! edge value and slope by computing the parabola as determined by
195  ! the right edge value and slope and the boundary cell average
196  u0_l = 3.0 * u0 + 0.5 * u1_r - 2.0 * u0_r
197 
198  ! Apply the traditional PPM limiter
199  exp1 = (u0_r - u0_l) * (u0 - 0.5*(u0_l+u0_r))
200  exp2 = (u0_r - u0_l) * (u0_r - u0_l) / 6.0
201 
202  if ( exp1 .GT. exp2 ) then
203  u0_l = 3.0 * u0 - 2.0 * u0_r
204  end if
205 
206  if ( exp1 .LT. -exp2 ) then
207  u0_r = 3.0 * u0 - 2.0 * u0_l
208  end if
209 
210  ppoly_e(i0,1) = u0_l
211  ppoly_e(i0,2) = u0_r
212 
213  a = u0_l
214  b = 6.0 * u0 - 4.0 * u0_l - 2.0 * u0_r
215  c = 3.0 * ( u0_r + u0_l - 2.0 * u0 )
216 
217  ppoly_coefficients(i0,1) = a
218  ppoly_coefficients(i0,2) = b
219  ppoly_coefficients(i0,3) = c
220 
221  ! ----- Right boundary -----
222  i0 = n-1
223  i1 = n
224  h0 = h(i0)
225  h1 = h(i1)
226  u0 = u(i0)
227  u1 = u(i1)
228 
229  ! Compute the right edge slope in neighboring cell and express it in
230  ! the global coordinate system
231  b = ppoly_coefficients(i0,2)
232  c = ppoly_coefficients(i0,3)
233  u1_l = (b + 2*c) ! derivative evaluated at xi = 1.0
234  u1_l = u1_l * ((h1+h_neglect)/(h0+h_neglect))
235 
236  ! Limit the left slope by the PLM limited slope
237  slope = 2.0 * ( u1 - u0 )
238  if ( abs(u1_l) .GT. abs(slope) ) then
239  u1_l = slope
240  end if
241 
242  ! The left edge value in the boundary cell is taken to be the right
243  ! edge value in the neighboring cell
244  u0_l = ppoly_e(i0,2)
245 
246  ! Given the left edge value and slope, we determine the right
247  ! edge value and slope by computing the parabola as determined by
248  ! the left edge value and slope and the boundary cell average
249  u0_r = 3.0 * u1 - 0.5 * u1_l - 2.0 * u0_l
250 
251  ! Apply the traditional PPM limiter
252  exp1 = (u0_r - u0_l) * (u1 - 0.5*(u0_l+u0_r))
253  exp2 = (u0_r - u0_l) * (u0_r - u0_l) / 6.0
254 
255  if ( exp1 .GT. exp2 ) then
256  u0_l = 3.0 * u1 - 2.0 * u0_r
257  end if
258 
259  if ( exp1 .LT. -exp2 ) then
260  u0_r = 3.0 * u1 - 2.0 * u0_l
261  end if
262 
263  ppoly_e(i1,1) = u0_l
264  ppoly_e(i1,2) = u0_r
265 
266  a = u0_l
267  b = 6.0 * u1 - 4.0 * u0_l - 2.0 * u0_r
268  c = 3.0 * ( u0_r + u0_l - 2.0 * u1 )
269 
270  ppoly_coefficients(i1,1) = a
271  ppoly_coefficients(i1,2) = b
272  ppoly_coefficients(i1,3) = c
273 
Here is the caller graph for this function:

◆ ppm_limiter_standard()

subroutine ppm_functions::ppm_limiter_standard ( integer, intent(in)  N,
real, dimension(n), intent(in)  h,
real, dimension(n), intent(in)  u,
real, dimension(n,2), intent(inout)  ppoly_E 
)
private

Adjusts edge values using the standard PPM limiter (Colella & Woodward, JCP 1984) after first checking that the edge values are bounded by neighbors cell averages and that the edge values are monotonic between cell averages.

Parameters
[in,out]ppoly_eEdge values

Definition at line 60 of file PPM_functions.F90.

References regrid_edge_values::bound_edge_values(), and regrid_edge_values::check_discontinuous_edge_values().

Referenced by ppm_reconstruction().

60  integer, intent(in) :: n ! Number of cells
61  real, dimension(N), intent(in) :: h ! Cell widths
62  real, dimension(N), intent(in) :: u ! Cell averages
63  real, dimension(N,2), intent(inout) :: ppoly_e !< Edge values
64  ! Local variables
65  integer :: k ! Loop index
66  real :: u_l, u_c, u_r ! Cell averages (left, center and right)
67  real :: edge_l, edge_r ! Edge values (left and right)
68  real :: expr1, expr2
69 
70  ! Bound edge values
71  call bound_edge_values( n, h, u, ppoly_e )
72 
73  ! Make discontinuous edge values monotonic
74  call check_discontinuous_edge_values( n, u, ppoly_e )
75 
76  ! Loop on interior cells to apply the standard
77  ! PPM limiter (Colella & Woodward, JCP 84)
78  do k = 2,n-1
79 
80  ! Get cell averages
81  u_l = u(k-1)
82  u_c = u(k)
83  u_r = u(k+1)
84 
85  edge_l = ppoly_e(k,1)
86  edge_r = ppoly_e(k,2)
87 
88  if ( (u_r - u_c)*(u_c - u_l) <= 0.0) then
89  ! Flatten extremum
90  edge_l = u_c
91  edge_r = u_c
92  else
93  expr1 = 3.0 * (edge_r - edge_l) * ( (u_c - edge_l) + (u_c - edge_r))
94  expr2 = (edge_r - edge_l) * (edge_r - edge_l)
95  if ( expr1 > expr2 ) then
96  ! Place extremum at right edge of cell by adjusting left edge value
97  edge_l = u_c + 2.0 * ( u_c - edge_r )
98  edge_l = max( min( edge_l, max(u_l, u_c) ), min(u_l, u_c) ) ! In case of round off
99  elseif ( expr1 < -expr2 ) then
100  ! Place extremum at left edge of cell by adjusting right edge value
101  edge_r = u_c + 2.0 * ( u_c - edge_l )
102  edge_r = max( min( edge_r, max(u_r, u_c) ), min(u_r, u_c) ) ! In case of round off
103  endif
104  endif
105  ! This checks that the difference in edge values is representable
106  ! and avoids overshoot problems due to round off.
107  if ( abs( edge_r - edge_l )<max(1.e-60,epsilon(u_c)*abs(u_c)) ) then
108  edge_l = u_c
109  edge_r = u_c
110  endif
111 
112  ppoly_e(k,1) = edge_l
113  ppoly_e(k,2) = edge_r
114 
115  enddo ! end loop on interior cells
116 
117  ! PCM within boundary cells
118  ppoly_e(1,:) = u(1)
119  ppoly_e(n,:) = u(n)
120 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ppm_reconstruction()

subroutine, public ppm_functions::ppm_reconstruction ( integer, intent(in)  N,
real, dimension(n), intent(in)  h,
real, dimension(n), intent(in)  u,
real, dimension(n,2), intent(inout)  ppoly_E,
real, dimension(n,3), intent(inout)  ppoly_coefficients 
)

Builds quadratic polynomials coefficients from cell mean and edge values.

Parameters
[in]nNumber of cells
[in]hCell widths
[in]uCell averages
[in,out]ppoly_eEdge values
[in,out]ppoly_coefficientsPolynomial coefficients

Definition at line 29 of file PPM_functions.F90.

References ppm_limiter_standard().

Referenced by mom_remapping::build_reconstructions_1d(), mom_ale::pressure_gradient_ppm(), regrid_interp::regridding_set_ppolys(), and mom_remapping::remapping_unit_tests().

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  real, dimension(N,3), intent(inout) :: ppoly_coefficients !< Polynomial coefficients
34  ! Local variables
35  integer :: k ! Loop index
36  real :: edge_l, edge_r ! Edge values (left and right)
37 
38  ! PPM limiter
39  call ppm_limiter_standard( n, h, u, ppoly_e )
40 
41  ! Loop over all cells
42  do k = 1,n
43 
44  edge_l = ppoly_e(k,1)
45  edge_r = ppoly_e(k,2)
46 
47  ! Store polynomial coefficients
48  ppoly_coefficients(k,1) = edge_l
49  ppoly_coefficients(k,2) = 4.0 * ( u(k) - edge_l ) + 2.0 * ( u(k) - edge_r )
50  ppoly_coefficients(k,3) = 3.0 * ( ( edge_r - u(k) ) + ( edge_l - u(k) ) )
51 
52  enddo
53 
Here is the call graph for this function:
Here is the caller graph for this function:

Variable Documentation

◆ h_neglect

real, parameter ppm_functions::h_neglect = 1.E-30
private

A tiny width that is so small that adding it to cell widths does not change the value due to a computational representation. It is used to avoid division by zero.

Note
This is a dimensional parameter and should really include a unit conversion.

Definition at line 23 of file PPM_functions.F90.

Referenced by ppm_boundary_extrapolation().

23 real, parameter :: h_neglect = 1.e-30