MOM6
Functions/Subroutines | Variables
plm_functions Module Reference

Detailed Description

Piecewise linear reconstruction functions.

Date of creation: 2008.06.06 L. White

This module contains routines that handle one-dimensionnal finite volume reconstruction using the piecewise linear method (PLM).

Functions/Subroutines

subroutine, public plm_reconstruction (N, h, u, ppoly_E, ppoly_coef, h_neglect)
 Reconstruction by linear polynomials within each cell. More...
 
subroutine, public plm_boundary_extrapolation (N, h, u, ppoly_E, ppoly_coef, h_neglect)
 Reconstruction by linear polynomials within boundary cells. More...
 

Variables

real, parameter hneglect_dflt = 1.E-30
 Default negligible cell thickness. More...
 

Function/Subroutine Documentation

◆ plm_boundary_extrapolation()

subroutine, public plm_functions::plm_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_coef,
real, intent(in), optional  h_neglect 
)

Reconstruction by linear polynomials within boundary cells.

The left and right edge values in the left and right boundary cells, respectively, are estimated using a linear extrapolation within the cells.

This extrapolation is EXACT when the underlying profile is linear.

It is assumed that the size of the array 'u' is equal to the number of cells defining 'grid' and 'ppoly'. No consistency check is performed here.

Parameters
[in]nNumber of cells
[in]hcell widths (size N)
[in]ucell averages (size N)
[in,out]ppoly_eedge values of piecewise polynomials, with the same units as u.
[in,out]ppoly_coefcoefficients of piecewise polynomials, mainly with the same units as u.
[in]h_neglectA negligibly small width for the purpose of cell reconstructions in the same units as h

Definition at line 206 of file PLM_functions.F90.

206  integer, intent(in) :: N !< Number of cells
207  real, dimension(:), intent(in) :: h !< cell widths (size N)
208  real, dimension(:), intent(in) :: u !< cell averages (size N)
209  real, dimension(:,:), intent(inout) :: ppoly_E !< edge values of piecewise polynomials,
210  !! with the same units as u.
211  real, dimension(:,:), intent(inout) :: ppoly_coef !< coefficients of piecewise polynomials, mainly
212  !! with the same units as u.
213  real, optional, intent(in) :: h_neglect !< A negligibly small width for
214  !! the purpose of cell reconstructions
215  !! in the same units as h
216 
217  ! Local variables
218  real :: u0, u1 ! cell averages
219  real :: h0, h1 ! corresponding cell widths
220  real :: slope ! retained PLM slope
221  real :: hNeglect
222 
223  hneglect = hneglect_dflt ; if (present(h_neglect)) hneglect = h_neglect
224 
225  ! -----------------------------------------
226  ! Left edge value in the left boundary cell
227  ! -----------------------------------------
228  h0 = h(1) + hneglect
229  h1 = h(2) + hneglect
230 
231  u0 = u(1)
232  u1 = u(2)
233 
234  ! The h2 scheme is used to compute the right edge value
235  ppoly_e(1,2) = (u0*h1 + u1*h0) / (h0 + h1)
236 
237  ! The standard PLM slope is computed as a first estimate for the
238  ! reconstruction within the cell
239  slope = 2.0 * ( ppoly_e(1,2) - u0 )
240 
241  ppoly_e(1,1) = u0 - 0.5 * slope
242  ppoly_e(1,2) = u0 + 0.5 * slope
243 
244  ppoly_coef(1,1) = ppoly_e(1,1)
245  ppoly_coef(1,2) = ppoly_e(1,2) - ppoly_e(1,1)
246 
247  ! ------------------------------------------
248  ! Right edge value in the left boundary cell
249  ! ------------------------------------------
250  h0 = h(n-1) + hneglect
251  h1 = h(n) + hneglect
252 
253  u0 = u(n-1)
254  u1 = u(n)
255 
256  ! The h2 scheme is used to compute the right edge value
257  ppoly_e(n,1) = (u0*h1 + u1*h0) / (h0 + h1)
258 
259  ! The standard PLM slope is computed as a first estimate for the
260  ! reconstruction within the cell
261  slope = 2.0 * ( u1 - ppoly_e(n,1) )
262 
263  ppoly_e(n,1) = u1 - 0.5 * slope
264  ppoly_e(n,2) = u1 + 0.5 * slope
265 
266  ppoly_coef(n,1) = ppoly_e(n,1)
267  ppoly_coef(n,2) = ppoly_e(n,2) - ppoly_e(n,1)
268 

References hneglect_dflt.

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

Here is the caller graph for this function:

◆ plm_reconstruction()

subroutine, public plm_functions::plm_reconstruction ( 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_coef,
real, intent(in), optional  h_neglect 
)

Reconstruction by linear polynomials within each cell.

It is assumed that the size of the array 'u' is equal to the number of cells defining 'grid' and 'ppoly'. No consistency check is performed here.

Parameters
[in]nNumber of cells
[in]hcell widths (size N)
[in]ucell averages (size N)
[in,out]ppoly_eedge values of piecewise polynomials, with the same units as u.
[in,out]ppoly_coefcoefficients of piecewise polynomials, mainly with the same units as u.
[in]h_neglectA negligibly small width for the purpose of cell reconstructions in the same units as h

Definition at line 19 of file PLM_functions.F90.

19  integer, intent(in) :: N !< Number of cells
20  real, dimension(:), intent(in) :: h !< cell widths (size N)
21  real, dimension(:), intent(in) :: u !< cell averages (size N)
22  real, dimension(:,:), intent(inout) :: ppoly_E !< edge values of piecewise polynomials,
23  !! with the same units as u.
24  real, dimension(:,:), intent(inout) :: ppoly_coef !< coefficients of piecewise polynomials, mainly
25  !! with the same units as u.
26  real, optional, intent(in) :: h_neglect !< A negligibly small width for
27  !! the purpose of cell reconstructions
28  !! in the same units as h
29 
30  ! Local variables
31  integer :: k ! loop index
32  real :: u_l, u_c, u_r ! left, center and right cell averages
33  real :: h_l, h_c, h_r, h_cn ! left, center and right cell widths
34  real :: sigma_l, sigma_c, sigma_r ! left, center and right
35  ! van Leer slopes
36  real :: slope ! retained PLM slope
37  real :: a, b ! auxiliary variables
38  real :: u_min, u_max, e_l, e_r, edge
39  real :: almost_one, almost_two
40  real, dimension(N) :: slp, mslp
41  real :: hNeglect
42 
43  hneglect = hneglect_dflt ; if (present(h_neglect)) hneglect = h_neglect
44 
45  almost_one = 1. - epsilon(slope)
46  almost_two = 2. * almost_one
47 
48  ! Loop on interior cells
49  do k = 2,n-1
50 
51  ! Get cell averages
52  u_l = u(k-1) ; u_c = u(k) ; u_r = u(k+1)
53 
54  ! Get cell widths
55  h_l = h(k-1) ; h_c = h(k) ; h_r = h(k+1)
56  h_cn = max( h_c, hneglect ) ! To avoid division by zero
57 
58  ! Side differences
59  sigma_r = u_r - u_c
60  sigma_l = u_c - u_l
61 
62  ! This is the second order slope given by equation 1.7 of
63  ! Piecewise Parabolic Method, Colella and Woodward (1984),
64  ! http://dx.doi.org/10.1016/0021-991(84)90143-8.
65  ! For uniform resolution it simplifies to ( u_r - u_l )/2 .
66  ! sigma_c = ( h_c / ( h_cn + ( h_l + h_r ) ) ) * ( &
67  ! ( 2.*h_l + h_c ) / ( h_r + h_cn ) * sigma_r &
68  ! + ( 2.*h_r + h_c ) / ( h_l + h_cn ) * sigma_l )
69 
70  ! This is the original estimate of the second order slope from Laurent
71  ! but multiplied by h_c
72  sigma_c = 2.0 * ( u_r - u_l ) * ( h_c / ( h_l + 2.0*h_c + h_r + hneglect) )
73 
74  if ( (sigma_l * sigma_r) > 0.0 ) then
75  ! This limits the slope so that the edge values are bounded by the
76  ! two cell averages spanning the edge.
77  u_min = min( u_l, u_c, u_r )
78  u_max = max( u_l, u_c, u_r )
79  slope = sign( min( abs(sigma_c), 2.*min( u_c - u_min, u_max - u_c ) ), sigma_c )
80  else
81  ! Extrema in the mean values require a PCM reconstruction avoid generating
82  ! larger extreme values.
83  slope = 0.0
84  endif
85 
86  ! This block tests to see if roundoff causes edge values to be out of bounds
87  u_min = min( u_l, u_c, u_r )
88  u_max = max( u_l, u_c, u_r )
89  if (u_c - 0.5*abs(slope) < u_min .or. u_c + 0.5*abs(slope) > u_max) then
90  slope = slope * almost_one
91  endif
92 
93  ! An attempt to avoid inconsistency when the values become unrepresentable.
94  if (abs(slope) < 1.e-140) slope = 0.
95 
96  ! Safety check - this block really should not be needed ...
97 ! if (u_c - 0.5*abs(slope) < u_min .or. u_c + 0.5*abs(slope) > u_max) then
98 ! write(0,*) 'l,c,r=',u_l,u_c,u_r
99 ! write(0,*) 'min,max=',u_min,u_max
100 ! write(0,*) 'slp=',slope
101 ! sigma_l = u_c-0.5*abs(slope)
102 ! sigma_r = u_c+0.5*abs(slope)
103 ! write(0,*) 'lo,hi=',sigma_l,sigma_r
104 ! write(0,*) 'elo,ehi=',sigma_l-u_min,sigma_r-u_max
105 ! stop 'Limiter failed!'
106 ! endif
107 
108  slp(k) = slope
109  ppoly_e(k,1) = u_c - 0.5 * slope
110  ppoly_e(k,2) = u_c + 0.5 * slope
111 
112  enddo ! end loop on interior cells
113 
114  ! Boundary cells use PCM. Extrapolation is handled in a later routine.
115  slp(1) = 0.
116  ppoly_e(1,2) = u(1)
117  slp(n) = 0.
118  ppoly_e(n,1) = u(n)
119 
120  ! This loop adjusts the slope so that edge values are monotonic.
121  do k = 2, n-1
122  u_l = u(k-1) ; u_c = u(k) ; u_r = u(k+1)
123  e_r = ppoly_e(k-1,2) ! Right edge from cell k-1
124  e_l = ppoly_e(k+1,1) ! Left edge from cell k
125  mslp(k) = abs(slp(k))
126  u_min = min(e_r, u_c)
127  u_max = max(e_r, u_c)
128  edge = u_c - 0.5 * slp(k)
129  if ( ( edge - e_r ) * ( u_c - edge ) < 0. ) then
130  edge = 0.5 * ( edge + e_r ) ! * almost_one?
131  mslp(k) = min( mslp(k), abs( edge - u_c ) * almost_two )
132  endif
133  edge = u_c + 0.5 * slp(k)
134  if ( ( edge - u_c ) * ( e_l - edge ) < 0. ) then
135  edge = 0.5 * ( edge + e_l ) ! * almost_one?
136  mslp(k) = min( mslp(k), abs( edge - u_c ) * almost_two )
137  endif
138  enddo ! end loop on interior cells
139  mslp(1) = 0.
140  mslp(n) = 0.
141 
142  ! Check that the above adjustment worked
143 ! do K = 2, N-1
144 ! u_r = u(k-1) + 0.5 * sign( mslp(k-1), slp(k-1) ) ! Right edge from cell k-1
145 ! u_l = u(k) - 0.5 * sign( mslp(k), slp(k) ) ! Left edge from cell k
146 ! if ( (u(k)-u(k-1)) * (u_l-u_r) < 0. ) then
147 ! stop 'Adjustment failed!'
148 ! endif
149 ! enddo ! end loop on interior cells
150 
151  ! Store and return edge values and polynomial coefficients.
152  ppoly_e(1,1) = u(1)
153  ppoly_e(1,2) = u(1)
154  ppoly_coef(1,1) = u(1)
155  ppoly_coef(1,2) = 0.
156  do k = 2, n-1
157  slope = sign( mslp(k), slp(k) )
158  u_l = u(k) - 0.5 * slope ! Left edge value of cell k
159  u_r = u(k) + 0.5 * slope ! Right edge value of cell k
160 
161  ! Check that final edge values are bounded
162  u_min = min( u(k-1), u(k) )
163  u_max = max( u(k-1), u(k) )
164  if (u_l<u_min .or. u_l>u_max) then
165  write(0,*) 'u(k-1)=',u(k-1),'u(k)=',u(k),'slp=',slp(k),'u_l=',u_l
166  stop 'Left edge out of bounds'
167  endif
168  u_min = min( u(k+1), u(k) )
169  u_max = max( u(k+1), u(k) )
170  if (u_r<u_min .or. u_r>u_max) then
171  write(0,*) 'u(k)=',u(k),'u(k+1)=',u(k+1),'slp=',slp(k),'u_r=',u_r
172  stop 'Right edge out of bounds'
173  endif
174 
175  ppoly_e(k,1) = u_l
176  ppoly_e(k,2) = u_r
177  ppoly_coef(k,1) = u_l
178  ppoly_coef(k,2) = ( u_r - u_l )
179  ! Check to see if this evaluation of the polynomial at x=1 would be
180  ! monotonic w.r.t. the next cell's edge value. If not, scale back!
181  edge = ppoly_coef(k,2) + ppoly_coef(k,1)
182  e_r = u(k+1) - 0.5 * sign( mslp(k+1), slp(k+1) )
183  if ( (edge-u(k))*(e_r-edge)<0.) then
184  ppoly_coef(k,2) = ppoly_coef(k,2) * almost_one
185  endif
186  enddo
187  ppoly_e(n,1) = u(n)
188  ppoly_e(n,2) = u(n)
189  ppoly_coef(n,1) = u(n)
190  ppoly_coef(n,2) = 0.
191 

References hneglect_dflt.

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

Here is the caller graph for this function:

Variable Documentation

◆ hneglect_dflt

real, parameter plm_functions::hneglect_dflt = 1.E-30
private

Default negligible cell thickness.

Definition at line 10 of file PLM_functions.F90.

10 real, parameter :: hNeglect_dflt = 1.e-30 !< Default negligible cell thickness

Referenced by plm_boundary_extrapolation(), and plm_reconstruction().