MOM6
regrid_edge_values.F90
1 !> Edge value estimation for high-order resconstruction
2 module regrid_edge_values
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
6 use regrid_solvers, only : solve_linear_system, solve_tridiagonal_system
7 use polynomial_functions, only : evaluation_polynomial
8 
9 implicit none ; private
10 
11 ! -----------------------------------------------------------------------------
12 ! The following routines are visible to the outside world
13 ! -----------------------------------------------------------------------------
14 public bound_edge_values
15 public average_discontinuous_edge_values
16 public check_discontinuous_edge_values
17 public edge_values_explicit_h2
18 public edge_values_explicit_h4
19 public edge_values_implicit_h4
20 public edge_values_implicit_h6
21 
22 #undef __DO_SAFETY_CHECKS__
23 
24 ! The following parameters are used to avoid singular matrices for boundary
25 ! extrapolation. The are needed only in the case where thicknesses vanish
26 ! to a small enough values such that the eigenvalues of the matrix can not
27 ! be separated.
28 ! Specifying a dimensional parameter value, as is done here, is a terrible idea.
29 real, parameter :: hneglect_edge_dflt = 1.e-10 !< The default value for cut-off minimum
30  !! thickness for sum(h) in edge value inversions
31 real, parameter :: hneglect_dflt = 1.e-30 !< The default value for cut-off minimum
32  !! thickness for sum(h) in other calculations
33 real, parameter :: hminfrac = 1.e-5 !< A minimum fraction for min(h)/sum(h)
34 
35 contains
36 
37 !> Bound edge values by neighboring cell averages
38 !!
39 !! In this routine, we loop on all cells to bound their left and right
40 !! edge values by the cell averages. That is, the left edge value must lie
41 !! between the left cell average and the central cell average. A similar
42 !! reasoning applies to the right edge values.
43 !!
44 !! Both boundary edge values are set equal to the boundary cell averages.
45 !! Any extrapolation scheme is applied after this routine has been called.
46 !! Therefore, boundary cells are treated as if they were local extrama.
47 subroutine bound_edge_values( N, h, u, edge_val, h_neglect )
48  integer, intent(in) :: n !< Number of cells
49  real, dimension(:), intent(in) :: h !< cell widths (size N)
50  real, dimension(:), intent(in) :: u !< cell average properties (size N)
51  real, dimension(:,:), intent(inout) :: edge_val !< Potentially modified edge values,
52  !! with the same units as u.
53  real, optional, intent(in) :: h_neglect !< A negligibly small width
54  !! in the same units as h.
55  ! Local variables
56  integer :: k ! loop index
57  integer :: k0, k1, k2
58  real :: h_l, h_c, h_r
59  real :: u_l, u_c, u_r
60  real :: u0_l, u0_r
61  real :: sigma_l, sigma_c, sigma_r ! left, center and right
62  ! van Leer slopes
63  real :: slope ! retained PLM slope
64 
65  real :: hneglect ! A negligible thicness in the same units as h.
66 
67  hneglect = hneglect_dflt ; if (present(h_neglect)) hneglect = h_neglect
68 
69  ! Loop on cells to bound edge value
70  do k = 1,n
71 
72  ! For the sake of bounding boundary edge values, the left neighbor
73  ! of the left boundary cell is assumed to be the same as the left
74  ! boundary cell and the right neighbor of the right boundary cell
75  ! is assumed to be the same as the right boundary cell. This
76  ! effectively makes boundary cells look like extrema.
77  if ( k == 1 ) then
78  k0 = 1
79  k1 = 1
80  k2 = 2
81  elseif ( k == n ) then
82  k0 = n-1
83  k1 = n
84  k2 = n
85  else
86  k0 = k-1
87  k1 = k
88  k2 = k+1
89  endif
90 
91  ! All cells can now be treated equally
92  h_l = h(k0)
93  h_c = h(k1)
94  h_r = h(k2)
95 
96  u_l = u(k0)
97  u_c = u(k1)
98  u_r = u(k2)
99 
100  u0_l = edge_val(k,1)
101  u0_r = edge_val(k,2)
102 
103  sigma_l = 2.0 * ( u_c - u_l ) / ( h_c + hneglect )
104  sigma_c = 2.0 * ( u_r - u_l ) / ( h_l + 2.0*h_c + h_r + hneglect )
105  sigma_r = 2.0 * ( u_r - u_c ) / ( h_c + hneglect )
106 
107  if ( (sigma_l * sigma_r) > 0.0 ) then
108  slope = sign( min(abs(sigma_l),abs(sigma_c),abs(sigma_r)), sigma_c )
109  else
110  slope = 0.0
111  endif
112 
113  ! The limiter must be used in the local coordinate system to each cell.
114  ! Hence, we must multiply the slope by h1. The multiplication by 0.5 is
115  ! simply a way to make it useable in the limiter (cfr White and Adcroft
116  ! JCP 2008 Eqs 19 and 20)
117  slope = slope * h_c * 0.5
118 
119  if ( (u_l-u0_l)*(u0_l-u_c) < 0.0 ) then
120  u0_l = u_c - sign( min( abs(slope), abs(u0_l-u_c) ), slope )
121  endif
122 
123  if ( (u_r-u0_r)*(u0_r-u_c) < 0.0 ) then
124  u0_r = u_c + sign( min( abs(slope), abs(u0_r-u_c) ), slope )
125  endif
126 
127  ! Finally bound by neighboring cell means in case of round off
128  u0_l = max( min( u0_l, max(u_l, u_c) ), min(u_l, u_c) )
129  u0_r = max( min( u0_r, max(u_r, u_c) ), min(u_r, u_c) )
130 
131  ! Store edge values
132  edge_val(k,1) = u0_l
133  edge_val(k,2) = u0_r
134 
135  enddo ! loop on interior edges
136 
137 end subroutine bound_edge_values
138 
139 !> Replace discontinuous collocated edge values with their average
140 !!
141 !! For each interior edge, check whether the edge values are discontinuous.
142 !! If so, compute the average and replace the edge values by the average.
143 subroutine average_discontinuous_edge_values( N, edge_val )
144  integer, intent(in) :: n !< Number of cells
145  real, dimension(:,:), intent(inout) :: edge_val !< Edge values that may be modified
146  !! the second index size is 2.
147  ! Local variables
148  integer :: k ! loop index
149  real :: u0_minus ! left value at given edge
150  real :: u0_plus ! right value at given edge
151  real :: u0_avg ! avg value at given edge
152 
153  ! Loop on interior edges
154  do k = 1,n-1
155 
156  ! Edge value on the left of the edge
157  u0_minus = edge_val(k,2)
158 
159  ! Edge value on the right of the edge
160  u0_plus = edge_val(k+1,1)
161 
162  if ( u0_minus /= u0_plus ) then
163  u0_avg = 0.5 * ( u0_minus + u0_plus )
164  edge_val(k,2) = u0_avg
165  edge_val(k+1,1) = u0_avg
166  endif
167 
168  enddo ! end loop on interior edges
169 
170 end subroutine average_discontinuous_edge_values
171 
172 !> Check discontinuous edge values and replace them with their average if not monotonic
173 !!
174 !! For each interior edge, check whether the edge values are discontinuous.
175 !! If so and if they are not monotonic, replace each edge value by their average.
176 subroutine check_discontinuous_edge_values( N, u, edge_val )
177  integer, intent(in) :: n !< Number of cells
178  real, dimension(:), intent(in) :: u !< cell averages (size N)
179  real, dimension(:,:), intent(inout) :: edge_val !< Cell edge values with the same units as u.
180  ! Local variables
181  integer :: k ! loop index
182  real :: u0_minus ! left value at given edge
183  real :: u0_plus ! right value at given edge
184  real :: um_minus ! left cell average
185  real :: um_plus ! right cell average
186  real :: u0_avg ! avg value at given edge
187 
188  ! Loop on interior cells
189  do k = 1,n-1
190 
191  ! Edge value on the left of the edge
192  u0_minus = edge_val(k,2)
193 
194  ! Edge value on the right of the edge
195  u0_plus = edge_val(k+1,1)
196 
197  ! Left cell average
198  um_minus = u(k)
199 
200  ! Right cell average
201  um_plus = u(k+1)
202 
203  if ( (u0_plus - u0_minus)*(um_plus - um_minus) < 0.0 ) then
204  u0_avg = 0.5 * ( u0_minus + u0_plus )
205  u0_avg = max( min( u0_avg, max(um_minus, um_plus) ), min(um_minus, um_plus) )
206  edge_val(k,2) = u0_avg
207  edge_val(k+1,1) = u0_avg
208  endif
209 
210  enddo ! end loop on interior edges
211 
212 end subroutine check_discontinuous_edge_values
213 
214 
215 !> Compute h2 edge values (explicit second order accurate)
216 !! in the same units as h.
217 !
218 !! Compute edge values based on second-order explicit estimates.
219 !! These estimates are based on a straight line spanning two cells and evaluated
220 !! at the location of the middle edge. An interpolant spanning cells
221 !! k-1 and k is evaluated at edge k-1/2. The estimate for each edge is unique.
222 !!
223 !! k-1 k
224 !! ..--o------o------o--..
225 !! k-1/2
226 !!
227 !! Boundary edge values are set to be equal to the boundary cell averages.
228 subroutine edge_values_explicit_h2( N, h, u, edge_val, h_neglect )
229  integer, intent(in) :: n !< Number of cells
230  real, dimension(:), intent(in) :: h !< cell widths (size N)
231  real, dimension(:), intent(in) :: u !< cell average properties (size N)
232  real, dimension(:,:), intent(inout) :: edge_val !< Returned edge values, with the
233  !! same units as u; the second index size is 2.
234  real, optional, intent(in) :: h_neglect !< A negligibly small width
235  ! Local variables
236  integer :: k ! loop index
237  real :: h0, h1 ! cell widths
238  real :: u0, u1 ! cell averages
239  real :: hneglect ! A negligible thicness in the same units as h.
240 
241  hneglect = hneglect_edge_dflt ; if (present(h_neglect)) hneglect = h_neglect
242 
243  ! Loop on interior cells
244  do k = 2,n
245 
246  h0 = h(k-1)
247  h1 = h(k)
248 
249  ! Avoid singularities when h0+h1=0
250  if (h0+h1==0.) then
251  h0 = hneglect
252  h1 = hneglect
253  endif
254 
255  u0 = u(k-1)
256  u1 = u(k)
257 
258  ! Compute left edge value
259  edge_val(k,1) = ( u0*h1 + u1*h0 ) / ( h0 + h1 )
260 
261  ! Left edge value of the current cell is equal to right edge
262  ! value of left cell
263  edge_val(k-1,2) = edge_val(k,1)
264 
265  enddo ! end loop on interior cells
266 
267  ! Boundary edge values are simply equal to the boundary cell averages
268  edge_val(1,1) = u(1)
269  edge_val(n,2) = u(n)
270 
271 end subroutine edge_values_explicit_h2
272 
273 !> Compute h4 edge values (explicit fourth order accurate)
274 !! in the same units as h.
275 !!
276 !! Compute edge values based on fourth-order explicit estimates.
277 !! These estimates are based on a cubic interpolant spanning four cells
278 !! and evaluated at the location of the middle edge. An interpolant spanning
279 !! cells i-2, i-1, i and i+1 is evaluated at edge i-1/2. The estimate for
280 !! each edge is unique.
281 !!
282 !! i-2 i-1 i i+1
283 !! ..--o------o------o------o------o--..
284 !! i-1/2
285 !!
286 !! The first two edge values are estimated by evaluating the first available
287 !! cubic interpolant, i.e., the interpolant spanning cells 1, 2, 3 and 4.
288 !! Similarly, the last two edge values are estimated by evaluating the last
289 !! available interpolant.
290 !!
291 !! For this fourth-order scheme, at least four cells must exist.
292 subroutine edge_values_explicit_h4( N, h, u, edge_val, h_neglect )
293  integer, intent(in) :: n !< Number of cells
294  real, dimension(:), intent(in) :: h !< cell widths (size N)
295  real, dimension(:), intent(in) :: u !< cell average properties (size N)
296  real, dimension(:,:), intent(inout) :: edge_val !< Returned edge values, with the
297  !! same units as u; the second index size is 2.
298  real, optional, intent(in) :: h_neglect !< A negligibly small width
299  ! Local variables
300  integer :: i, j
301  real :: u0, u1, u2, u3
302  real :: h0, h1, h2, h3
303  real :: f1, f2, f3 ! auxiliary variables
304  real :: e ! edge value
305  real, dimension(5) :: x ! used to compute edge
306  real, dimension(4,4) :: a ! values near the boundaries
307  real, dimension(4) :: b, c
308  real :: hneglect ! A negligible thicness in the same units as h.
309 
310  hneglect = hneglect_edge_dflt ; if (present(h_neglect)) hneglect = h_neglect
311 
312  ! Loop on interior cells
313  do i = 3,n-1
314 
315  h0 = h(i-2)
316  h1 = h(i-1)
317  h2 = h(i)
318  h3 = h(i+1)
319 
320  ! Avoid singularities when consecutive pairs of h vanish
321  if (h0+h1==0. .or. h1+h2==0. .or. h2+h3==0.) then
322  f1 = max( hneglect, h0+h1+h2+h3 )
323  h0 = max( hminfrac*f1, h(i-2) )
324  h1 = max( hminfrac*f1, h(i-1) )
325  h2 = max( hminfrac*f1, h(i) )
326  h3 = max( hminfrac*f1, h(i+1) )
327  endif
328 
329  u0 = u(i-2)
330  u1 = u(i-1)
331  u2 = u(i)
332  u3 = u(i+1)
333 
334  f1 = (h0+h1) * (h2+h3) / (h1+h2)
335  f2 = u1 * h2 + u2 * h1
336  f3 = 1.0 / (h0+h1+h2) + 1.0 / (h1+h2+h3)
337 
338  e = f1 * f2 * f3
339 
340  f1 = h2 * (h2+h3) / ( (h0+h1+h2)*(h0+h1) )
341  f2 = u1*(h0+2.0*h1) - u0*h1
342 
343  e = e + f1*f2
344 
345  f1 = h1 * (h0+h1) / ( (h1+h2+h3)*(h2+h3) )
346  f2 = u2*(2.0*h2+h3) - u3*h2
347 
348  e = e + f1*f2
349 
350  e = e / ( h0 + h1 + h2 + h3)
351 
352  edge_val(i,1) = e
353  edge_val(i-1,2) = e
354 
355 #ifdef __DO_SAFETY_CHECKS__
356  if (e /= e) then
357  write(0,*) 'NaN in explicit_edge_h4 at k=',i
358  write(0,*) 'u0-u3=',u0,u1,u2,u3
359  write(0,*) 'h0-h3=',h0,h1,h2,h3
360  write(0,*) 'f1-f3=',f1,f2,f3
361  stop 'Nan during edge_values_explicit_h4'
362  endif
363 #endif
364 
365  enddo ! end loop on interior cells
366 
367  ! Determine first two edge values
368  f1 = max( hneglect, hminfrac*sum(h(1:4)) )
369  x(1) = 0.0
370  do i = 2,5
371  x(i) = x(i-1) + max(f1, h(i-1))
372  enddo
373 
374  do i = 1,4
375 
376  do j = 1,4
377  a(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / real(j)
378  enddo
379 
380  b(i) = u(i) * max(f1, h(i) )
381 
382  enddo
383 
384  call solve_linear_system( a, b, c, 4 )
385 
386  ! First edge value
387  edge_val(1,1) = evaluation_polynomial( c, 4, x(1) )
388 
389  ! Second edge value
390  edge_val(1,2) = evaluation_polynomial( c, 4, x(2) )
391  edge_val(2,1) = edge_val(1,2)
392 
393 #ifdef __DO_SAFETY_CHECKS__
394  if (edge_val(1,1) /= edge_val(1,1) .or. edge_val(1,2) /= edge_val(1,2)) then
395  write(0,*) 'NaN in explicit_edge_h4 at k=',1
396  write(0,*) 'A=',a
397  write(0,*) 'B=',b
398  write(0,*) 'C=',c
399  write(0,*) 'h(1:4)=',h(1:4)
400  write(0,*) 'x=',x
401  stop 'Nan during edge_values_explicit_h4'
402  endif
403 #endif
404 
405  ! Determine last two edge values
406  f1 = max( hneglect, hminfrac*sum(h(n-3:n)) )
407  x(1) = 0.0
408  do i = 2,5
409  x(i) = x(i-1) + max(f1, h(n-5+i))
410  enddo
411 
412  do i = 1,4
413 
414  do j = 1,4
415  a(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / real(j)
416  enddo
417 
418  b(i) = u(n-4+i) * max(f1, h(n-4+i) )
419 
420  enddo
421 
422  call solve_linear_system( a, b, c, 4 )
423 
424  ! Last edge value
425  edge_val(n,2) = evaluation_polynomial( c, 4, x(5) )
426 
427  ! Second to last edge value
428  edge_val(n,1) = evaluation_polynomial( c, 4, x(4) )
429  edge_val(n-1,2) = edge_val(n,1)
430 
431 #ifdef __DO_SAFETY_CHECKS__
432  if (edge_val(n,1) /= edge_val(n,1) .or. edge_val(n,2) /= edge_val(n,2)) then
433  write(0,*) 'NaN in explicit_edge_h4 at k=',n
434  write(0,*) 'A='
435  do i = 1,4
436  do j = 1,4
437  a(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / real(j)
438  enddo
439  write(0,*) a(i,:)
440  b(i) = u(n-4+i) * ( h(n-4+i) )
441  enddo
442  write(0,*) 'B=',b
443  write(0,*) 'C=',c
444  write(0,*) 'h(:N)=',h(n-3:n)
445  write(0,*) 'x=',x
446  stop 'Nan during edge_values_explicit_h4'
447  endif
448 #endif
449 
450 end subroutine edge_values_explicit_h4
451 
452 !> Compute ih4 edge values (implicit fourth order accurate)
453 !! in the same units as h.
454 !!
455 !! Compute edge values based on fourth-order implicit estimates.
456 !!
457 !! Fourth-order implicit estimates of edge values are based on a two-cell
458 !! stencil. A tridiagonal system is set up and is based on expressing the
459 !! edge values in terms of neighboring cell averages. The generic
460 !! relationship is
461 !!
462 !! \f[
463 !! \alpha u_{i-1/2} + u_{i+1/2} + \beta u_{i+3/2} = a \bar{u}_i + b \bar{u}_{i+1}
464 !! \f]
465 !!
466 !! and the stencil looks like this
467 !!
468 !! i i+1
469 !! ..--o------o------o--..
470 !! i-1/2 i+1/2 i+3/2
471 !!
472 !! In this routine, the coefficients \f$\alpha\f$, \f$\beta\f$, \f$a\f$ and \f$b\f$ are
473 !! computed, the tridiagonal system is built, boundary conditions are prescribed and
474 !! the system is solved to yield edge-value estimates.
475 !!
476 !! There are N+1 unknowns and we are able to write N-1 equations. The
477 !! boundary conditions close the system.
478 subroutine edge_values_implicit_h4( N, h, u, edge_val, h_neglect )
479  integer, intent(in) :: n !< Number of cells
480  real, dimension(:), intent(in) :: h !< cell widths (size N)
481  real, dimension(:), intent(in) :: u !< cell average properties (size N)
482  real, dimension(:,:), intent(inout) :: edge_val !< Returned edge values, with the
483  !! same units as u; the second index size is 2.
484  real, optional, intent(in) :: h_neglect !< A negligibly small width
485  ! Local variables
486  integer :: i, j ! loop indexes
487  real :: h0, h1 ! cell widths
488  real :: h0_2, h1_2, h0h1
489  real :: d2, d4
490  real :: alpha, beta ! stencil coefficients
491  real :: a, b
492  real, dimension(5) :: x ! system used to enforce
493  real, dimension(4,4) :: asys ! boundary conditions
494  real, dimension(4) :: bsys, csys
495  real, dimension(N+1) :: tri_l, & ! trid. system (lower diagonal)
496  tri_d, & ! trid. system (middle diagonal)
497  tri_u, & ! trid. system (upper diagonal)
498  tri_b, & ! trid. system (unknowns vector)
499  tri_x ! trid. system (rhs)
500  real :: hneglect ! A negligible thicness in the same units as h.
501 
502  hneglect = hneglect_edge_dflt ; if (present(h_neglect)) hneglect = h_neglect
503 
504  ! Loop on cells (except last one)
505  do i = 1,n-1
506 
507  ! Get cell widths
508  h0 = h(i)
509  h1 = h(i+1)
510 
511  ! Avoid singularities when h0+h1=0
512  if (h0+h1==0.) then
513  h0 = hneglect
514  h1 = hneglect
515  endif
516 
517  ! Auxiliary calculations
518  d2 = (h0 + h1) ** 2
519  d4 = d2 ** 2
520  h0h1 = h0 * h1
521  h0_2 = h0 * h0
522  h1_2 = h1 * h1
523 
524  ! Coefficients
525  alpha = h1_2 / d2
526  beta = h0_2 / d2
527  a = 2.0 * h1_2 * ( h1_2 + 2.0 * h0_2 + 3.0 * h0h1 ) / d4
528  b = 2.0 * h0_2 * ( h0_2 + 2.0 * h1_2 + 3.0 * h0h1 ) / d4
529 
530  tri_l(i+1) = alpha
531  tri_d(i+1) = 1.0
532  tri_u(i+1) = beta
533 
534  tri_b(i+1) = a * u(i) + b * u(i+1)
535 
536  enddo ! end loop on cells
537 
538  ! Boundary conditions: left boundary
539  h0 = max( hneglect, hminfrac*sum(h(1:4)) )
540  x(1) = 0.0
541  do i = 2,5
542  x(i) = x(i-1) + max( h0, h(i-1) )
543  enddo
544 
545  do i = 1,4
546 
547  do j = 1,4
548  asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j
549  enddo
550 
551  bsys(i) = u(i) * max( h0, h(i) )
552 
553  enddo
554 
555  call solve_linear_system( asys, bsys, csys, 4 )
556 
557  tri_d(1) = 1.0
558  tri_u(1) = 0.0
559  tri_b(1) = evaluation_polynomial( csys, 4, x(1) ) ! first edge value
560 
561  ! Boundary conditions: right boundary
562  h0 = max( hneglect, hminfrac*sum(h(n-3:n)) )
563  x(1) = 0.0
564  do i = 2,5
565  x(i) = x(i-1) + max( h0, h(n-5+i) )
566  enddo
567 
568  do i = 1,4
569 
570  do j = 1,4
571  asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j
572  enddo
573 
574  bsys(i) = u(n-4+i) * max( h0, h(n-4+i) )
575 
576  enddo
577 
578  call solve_linear_system( asys, bsys, csys, 4 )
579 
580  tri_l(n+1) = 0.0
581  tri_d(n+1) = 1.0
582  tri_b(n+1) = evaluation_polynomial( csys, 4, x(5) ) ! last edge value
583 
584  ! Solve tridiagonal system and assign edge values
585  call solve_tridiagonal_system( tri_l, tri_d, tri_u, tri_b, tri_x, n+1 )
586 
587  do i = 2,n
588  edge_val(i,1) = tri_x(i)
589  edge_val(i-1,2) = tri_x(i)
590  enddo
591  edge_val(1,1) = tri_x(1)
592  edge_val(n,2) = tri_x(n+1)
593 
594 end subroutine edge_values_implicit_h4
595 
596 !> Compute ih6 edge values (implicit sixth order accurate)
597  !! in the same units as h.
598 !!
599 !! Sixth-order implicit estimates of edge values are based on a four-cell,
600 !! three-edge stencil. A tridiagonal system is set up and is based on
601 !! expressing the edge values in terms of neighboring cell averages.
602 !!
603 !! The generic relationship is
604 !!
605 !! \f[
606 !! \alpha u_{i-1/2} + u_{i+1/2} + \beta u_{i+3/2} =
607 !! a \bar{u}_{i-1} + b \bar{u}_i + c \bar{u}_{i+1} + d \bar{u}_{i+2}
608 !! \f]
609 !!
610 !! and the stencil looks like this
611 !!
612 !! i-1 i i+1 i+2
613 !! ..--o------o------o------o------o--..
614 !! i-1/2 i+1/2 i+3/2
615 !!
616 !! In this routine, the coefficients \f$\alpha\f$, \f$\beta\f$, a, b, c and d are
617 !! computed, the tridiagonal system is built, boundary conditions are
618 !! prescribed and the system is solved to yield edge-value estimates.
619 !!
620 !! Note that the centered stencil only applies to edges 3 to N-1 (edges are
621 !! numbered 1 to n+1), which yields N-3 equations for N+1 unknowns. Two other
622 !! equations are written by using a right-biased stencil for edge 2 and a
623 !! left-biased stencil for edge N. The prescription of boundary conditions
624 !! (using sixth-order polynomials) closes the system.
625 !!
626 !! CAUTION: For each edge, in order to determine the coefficients of the
627 !! implicit expression, a 6x6 linear system is solved. This may
628 !! become computationally expensive if regridding is carried out
629 !! often. Figuring out closed-form expressions for these coefficients
630 !! on nonuniform meshes turned out to be intractable.
631 subroutine edge_values_implicit_h6( N, h, u, edge_val, h_neglect )
632  integer, intent(in) :: n !< Number of cells
633  real, dimension(:), intent(in) :: h !< cell widths (size N)
634  real, dimension(:), intent(in) :: u !< cell average properties (size N)
635  real, dimension(:,:), intent(inout) :: edge_val !< Returned edge values, with the
636  !! same units as u; the second index size is 2.
637  real, optional, intent(in) :: h_neglect !< A negligibly small width
638  ! Local variables
639  integer :: i, j, k ! loop indexes
640  real :: h0, h1, h2, h3 ! cell widths
641  real :: g, g_2, g_3 ! the following are
642  real :: g_4, g_5, g_6 ! auxiliary variables
643  real :: d2, d3, d4, d5, d6 ! to set up the systems
644  real :: n2, n3, n4, n5, n6 ! used to compute the
645  real :: h1_2, h2_2 ! the coefficients of the
646  real :: h1_3, h2_3 ! tridiagonal system
647  real :: h1_4, h2_4 ! ...
648  real :: h1_5, h2_5 ! ...
649  real :: h1_6, h2_6 ! ...
650  real :: h0ph1, h0ph1_2 ! ...
651  real :: h0ph1_3, h0ph1_4 ! ...
652  real :: h2ph3, h2ph3_2 ! ...
653  real :: h2ph3_3, h2ph3_4 ! ...
654  real :: h0ph1_5, h2ph3_5 ! ...
655  real :: alpha, beta ! stencil coefficients
656  real :: a, b, c, d ! "
657  real, dimension(7) :: x ! system used to enforce
658  real, dimension(6,6) :: asys ! boundary conditions
659  real, dimension(6) :: bsys, csys ! ...
660  real, dimension(N+1) :: tri_l, & ! trid. system (lower diagonal)
661  tri_d, & ! trid. system (middle diagonal)
662  tri_u, & ! trid. system (upper diagonal)
663  tri_b, & ! trid. system (unknowns vector)
664  tri_x ! trid. system (rhs)
665  real :: hneglect ! A negligible thicness in the same units as h.
666 
667  hneglect = hneglect_edge_dflt ; if (present(h_neglect)) hneglect = h_neglect
668 
669  ! Loop on cells (except last one)
670  do k = 2,n-2
671 
672  ! Cell widths
673  h0 = h(k-1)
674  h1 = h(k+0)
675  h2 = h(k+1)
676  h3 = h(k+2)
677 
678  ! Avoid singularities when h0=0 or h3=0
679  if (h0*h3==0.) then
680  g = max( hneglect, h0+h1+h2+h3 )
681  h0 = max( hminfrac*g, h0 )
682  h1 = max( hminfrac*g, h1 )
683  h2 = max( hminfrac*g, h2 )
684  h3 = max( hminfrac*g, h3 )
685  endif
686 
687  ! Auxiliary calculations
688  h1_2 = h1 * h1
689  h1_3 = h1_2 * h1
690  h1_4 = h1_2 * h1_2
691  h1_5 = h1_3 * h1_2
692  h1_6 = h1_3 * h1_3
693 
694  h2_2 = h2 * h2
695  h2_3 = h2_2 * h2
696  h2_4 = h2_2 * h2_2
697  h2_5 = h2_3 * h2_2
698  h2_6 = h2_3 * h2_3
699 
700  g = h0 + h1
701  g_2 = g * g
702  g_3 = g * g_2
703  g_4 = g_2 * g_2
704  g_5 = g_4 * g
705  g_6 = g_3 * g_3
706 
707  d2 = ( h1_2 - g_2 ) / h0
708  d3 = ( h1_3 - g_3 ) / h0
709  d4 = ( h1_4 - g_4 ) / h0
710  d5 = ( h1_5 - g_5 ) / h0
711  d6 = ( h1_6 - g_6 ) / h0
712 
713  g = h2 + h3
714  g_2 = g * g
715  g_3 = g * g_2
716  g_4 = g_2 * g_2
717  g_5 = g_4 * g
718  g_6 = g_3 * g_3
719 
720  n2 = ( g_2 - h2_2 ) / h3
721  n3 = ( g_3 - h2_3 ) / h3
722  n4 = ( g_4 - h2_4 ) / h3
723  n5 = ( g_5 - h2_5 ) / h3
724  n6 = ( g_6 - h2_6 ) / h3
725 
726  ! Compute matrix entries
727  asys(1,1) = 1.0
728  asys(1,2) = 1.0
729  asys(1,3) = -1.0
730  asys(1,4) = -1.0
731  asys(1,5) = -1.0
732  asys(1,6) = -1.0
733 
734  asys(2,1) = - h1
735  asys(2,2) = h2
736  asys(2,3) = -0.5 * d2
737  asys(2,4) = 0.5 * h1
738  asys(2,5) = -0.5 * h2
739  asys(2,6) = -0.5 * n2
740 
741  asys(3,1) = 0.5 * h1_2
742  asys(3,2) = 0.5 * h2_2
743  asys(3,3) = d3 / 6.0
744  asys(3,4) = - h1_2 / 6.0
745  asys(3,5) = - h2_2 / 6.0
746  asys(3,6) = - n3 / 6.0
747 
748  asys(4,1) = - h1_3 / 6.0
749  asys(4,2) = h2_3 / 6.0
750  asys(4,3) = - d4 / 24.0
751  asys(4,4) = h1_3 / 24.0
752  asys(4,5) = - h2_3 / 24.0
753  asys(4,6) = - n4 / 24.0
754 
755  asys(5,1) = h1_4 / 24.0
756  asys(5,2) = h2_4 / 24.0
757  asys(5,3) = d5 / 120.0
758  asys(5,4) = - h1_4 / 120.0
759  asys(5,5) = - h2_4 / 120.0
760  asys(5,6) = - n5 / 120.0
761 
762  asys(6,1) = - h1_5 / 120.0
763  asys(6,2) = h2_5 / 120.0
764  asys(6,3) = - d6 / 720.0
765  asys(6,4) = h1_5 / 720.0
766  asys(6,5) = - h2_5 / 720.0
767  asys(6,6) = - n6 / 720.0
768 
769  bsys(:) = (/ -1.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)
770 
771  call solve_linear_system( asys, bsys, csys, 6 )
772 
773  alpha = csys(1)
774  beta = csys(2)
775  a = csys(3)
776  b = csys(4)
777  c = csys(5)
778  d = csys(6)
779 
780  tri_l(k+1) = alpha
781  tri_d(k+1) = 1.0
782  tri_u(k+1) = beta
783  tri_b(k+1) = a * u(k-1) + b * u(k) + c * u(k+1) + d * u(k+2)
784 
785  enddo ! end loop on cells
786 
787  ! Use a right-biased stencil for the second row
788 
789  ! Cell widths
790  h0 = h(1)
791  h1 = h(2)
792  h2 = h(3)
793  h3 = h(4)
794 
795  ! Avoid singularities when h0=0 or h3=0
796  if (h0*h3==0.) then
797  g = max( hneglect, h0+h1+h2+h3 )
798  h0 = max( hminfrac*g, h0 )
799  h1 = max( hminfrac*g, h1 )
800  h2 = max( hminfrac*g, h2 )
801  h3 = max( hminfrac*g, h3 )
802  endif
803 
804  ! Auxiliary calculations
805  h1_2 = h1 * h1
806  h1_3 = h1_2 * h1
807  h1_4 = h1_2 * h1_2
808  h1_5 = h1_3 * h1_2
809  h1_6 = h1_3 * h1_3
810 
811  h2_2 = h2 * h2
812  h2_3 = h2_2 * h2
813  h2_4 = h2_2 * h2_2
814  h2_5 = h2_3 * h2_2
815  h2_6 = h2_3 * h2_3
816 
817  g = h0 + h1
818  g_2 = g * g
819  g_3 = g * g_2
820  g_4 = g_2 * g_2
821  g_5 = g_4 * g
822  g_6 = g_3 * g_3
823 
824  h0ph1 = h0 + h1
825  h0ph1_2 = h0ph1 * h0ph1
826  h0ph1_3 = h0ph1_2 * h0ph1
827  h0ph1_4 = h0ph1_2 * h0ph1_2
828  h0ph1_5 = h0ph1_3 * h0ph1_2
829 
830  d2 = ( h1_2 - g_2 ) / h0
831  d3 = ( h1_3 - g_3 ) / h0
832  d4 = ( h1_4 - g_4 ) / h0
833  d5 = ( h1_5 - g_5 ) / h0
834  d6 = ( h1_6 - g_6 ) / h0
835 
836  g = h2 + h3
837  g_2 = g * g
838  g_3 = g * g_2
839  g_4 = g_2 * g_2
840  g_5 = g_4 * g
841  g_6 = g_3 * g_3
842 
843  n2 = ( g_2 - h2_2 ) / h3
844  n3 = ( g_3 - h2_3 ) / h3
845  n4 = ( g_4 - h2_4 ) / h3
846  n5 = ( g_5 - h2_5 ) / h3
847  n6 = ( g_6 - h2_6 ) / h3
848 
849  ! Compute matrix entries
850  asys(1,1) = 1.0
851  asys(1,2) = 1.0
852  asys(1,3) = -1.0
853  asys(1,4) = -1.0
854  asys(1,5) = -1.0
855  asys(1,6) = -1.0
856 
857  asys(2,1) = - h0ph1
858  asys(2,2) = 0.0
859  asys(2,3) = -0.5 * d2
860  asys(2,4) = 0.5 * h1
861  asys(2,5) = -0.5 * h2
862  asys(2,6) = -0.5 * n2
863 
864  asys(3,1) = 0.5 * h0ph1_2
865  asys(3,2) = 0.0
866  asys(3,3) = d3 / 6.0
867  asys(3,4) = - h1_2 / 6.0
868  asys(3,5) = - h2_2 / 6.0
869  asys(3,6) = - n3 / 6.0
870 
871  asys(4,1) = - h0ph1_3 / 6.0
872  asys(4,2) = 0.0
873  asys(4,3) = - d4 / 24.0
874  asys(4,4) = h1_3 / 24.0
875  asys(4,5) = - h2_3 / 24.0
876  asys(4,6) = - n4 / 24.0
877 
878  asys(5,1) = h0ph1_4 / 24.0
879  asys(5,2) = 0.0
880  asys(5,3) = d5 / 120.0
881  asys(5,4) = - h1_4 / 120.0
882  asys(5,5) = - h2_4 / 120.0
883  asys(5,6) = - n5 / 120.0
884 
885  asys(6,1) = - h0ph1_5 / 120.0
886  asys(6,2) = 0.0
887  asys(6,3) = - d6 / 720.0
888  asys(6,4) = h1_5 / 720.0
889  asys(6,5) = - h2_5 / 720.0
890  asys(6,6) = - n6 / 720.0
891 
892  bsys(:) = (/ -1.0, h1, -0.5*h1_2, h1_3/6.0, -h1_4/24.0, h1_5/120.0 /)
893 
894  call solve_linear_system( asys, bsys, csys, 6 )
895 
896  alpha = csys(1)
897  beta = csys(2)
898  a = csys(3)
899  b = csys(4)
900  c = csys(5)
901  d = csys(6)
902 
903  tri_l(2) = alpha
904  tri_d(2) = 1.0
905  tri_u(2) = beta
906  tri_b(2) = a * u(1) + b * u(2) + c * u(3) + d * u(4)
907 
908  ! Boundary conditions: left boundary
909  g = max( hneglect, hminfrac*sum(h(1:6)) )
910  x(1) = 0.0
911  do i = 2,7
912  x(i) = x(i-1) + max( g, h(i-1) )
913  enddo
914 
915  do i = 1,6
916 
917  do j = 1,6
918  asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j
919  enddo
920 
921  bsys(i) = u(i) * max( g, h(i) )
922 
923  enddo
924 
925  call solve_linear_system( asys, bsys, csys, 6 )
926 
927  tri_l(1) = 0.0
928  tri_d(1) = 1.0
929  tri_u(1) = 0.0
930  tri_b(1) = evaluation_polynomial( csys, 6, x(1) ) ! first edge value
931 
932  ! Use a left-biased stencil for the second to last row
933 
934  ! Cell widths
935  h0 = h(n-3)
936  h1 = h(n-2)
937  h2 = h(n-1)
938  h3 = h(n)
939 
940  ! Avoid singularities when h0=0 or h3=0
941  if (h0*h3==0.) then
942  g = max( hneglect, h0+h1+h2+h3 )
943  h0 = max( hminfrac*g, h0 )
944  h1 = max( hminfrac*g, h1 )
945  h2 = max( hminfrac*g, h2 )
946  h3 = max( hminfrac*g, h3 )
947  endif
948 
949  ! Auxiliary calculations
950  h1_2 = h1 * h1
951  h1_3 = h1_2 * h1
952  h1_4 = h1_2 * h1_2
953  h1_5 = h1_3 * h1_2
954  h1_6 = h1_3 * h1_3
955 
956  h2_2 = h2 * h2
957  h2_3 = h2_2 * h2
958  h2_4 = h2_2 * h2_2
959  h2_5 = h2_3 * h2_2
960  h2_6 = h2_3 * h2_3
961 
962  g = h0 + h1
963  g_2 = g * g
964  g_3 = g * g_2
965  g_4 = g_2 * g_2
966  g_5 = g_4 * g
967  g_6 = g_3 * g_3
968 
969  h2ph3 = h2 + h3
970  h2ph3_2 = h2ph3 * h2ph3
971  h2ph3_3 = h2ph3_2 * h2ph3
972  h2ph3_4 = h2ph3_2 * h2ph3_2
973  h2ph3_5 = h2ph3_3 * h2ph3_2
974 
975  d2 = ( h1_2 - g_2 ) / h0
976  d3 = ( h1_3 - g_3 ) / h0
977  d4 = ( h1_4 - g_4 ) / h0
978  d5 = ( h1_5 - g_5 ) / h0
979  d6 = ( h1_6 - g_6 ) / h0
980 
981  g = h2 + h3
982  g_2 = g * g
983  g_3 = g * g_2
984  g_4 = g_2 * g_2
985  g_5 = g_4 * g
986  g_6 = g_3 * g_3
987 
988  n2 = ( g_2 - h2_2 ) / h3
989  n3 = ( g_3 - h2_3 ) / h3
990  n4 = ( g_4 - h2_4 ) / h3
991  n5 = ( g_5 - h2_5 ) / h3
992  n6 = ( g_6 - h2_6 ) / h3
993 
994  ! Compute matrix entries
995  asys(1,1) = 1.0
996  asys(1,2) = 1.0
997  asys(1,3) = -1.0
998  asys(1,4) = -1.0
999  asys(1,5) = -1.0
1000  asys(1,6) = -1.0
1001 
1002  asys(2,1) = 0.0
1003  asys(2,2) = h2ph3
1004  asys(2,3) = -0.5 * d2
1005  asys(2,4) = 0.5 * h1
1006  asys(2,5) = -0.5 * h2
1007  asys(2,6) = -0.5 * n2
1008 
1009  asys(3,1) = 0.0
1010  asys(3,2) = 0.5 * h2ph3_2
1011  asys(3,3) = d3 / 6.0
1012  asys(3,4) = - h1_2 / 6.0
1013  asys(3,5) = - h2_2 / 6.0
1014  asys(3,6) = - n3 / 6.0
1015 
1016  asys(4,1) = 0.0
1017  asys(4,2) = h2ph3_3 / 6.0
1018  asys(4,3) = - d4 / 24.0
1019  asys(4,4) = h1_3 / 24.0
1020  asys(4,5) = - h2_3 / 24.0
1021  asys(4,6) = - n4 / 24.0
1022 
1023  asys(5,1) = 0.0
1024  asys(5,2) = h2ph3_4 / 24.0
1025  asys(5,3) = d5 / 120.0
1026  asys(5,4) = - h1_4 / 120.0
1027  asys(5,5) = - h2_4 / 120.0
1028  asys(5,6) = - n5 / 120.0
1029 
1030  asys(6,1) = 0.0
1031  asys(6,2) = h2ph3_5 / 120.0
1032  asys(6,3) = - d6 / 720.0
1033  asys(6,4) = h1_5 / 720.0
1034  asys(6,5) = - h2_5 / 720.0
1035  asys(6,6) = - n6 / 720.0
1036 
1037  bsys(:) = (/ -1.0, -h2, -0.5*h2_2, -h2_3/6.0, -h2_4/24.0, -h2_5/120.0 /)
1038 
1039  call solve_linear_system( asys, bsys, csys, 6 )
1040 
1041  alpha = csys(1)
1042  beta = csys(2)
1043  a = csys(3)
1044  b = csys(4)
1045  c = csys(5)
1046  d = csys(6)
1047 
1048  tri_l(n) = alpha
1049  tri_d(n) = 1.0
1050  tri_u(n) = beta
1051  tri_b(n) = a * u(n-3) + b * u(n-2) + c * u(n-1) + d * u(n)
1052 
1053  ! Boundary conditions: right boundary
1054  g = max( hneglect, hminfrac*sum(h(n-5:n)) )
1055  x(1) = 0.0
1056  do i = 2,7
1057  x(i) = x(i-1) + max( g, h(n-7+i) )
1058  enddo
1059 
1060  do i = 1,6
1061 
1062  do j = 1,6
1063  asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j
1064  enddo
1065 
1066  bsys(i) = u(n-6+i) * max( g, h(n-6+i) )
1067 
1068  enddo
1069 
1070  call solve_linear_system( asys, bsys, csys, 6 )
1071 
1072  tri_l(n+1) = 0.0
1073  tri_d(n+1) = 1.0
1074  tri_u(n+1) = 0.0
1075  tri_b(n+1) = evaluation_polynomial( csys, 6, x(7) ) ! last edge value
1076 
1077  ! Solve tridiagonal system and assign edge values
1078  call solve_tridiagonal_system( tri_l, tri_d, tri_u, tri_b, tri_x, n+1 )
1079 
1080  do i = 2,n
1081  edge_val(i,1) = tri_x(i)
1082  edge_val(i-1,2) = tri_x(i)
1083  enddo
1084  edge_val(1,1) = tri_x(1)
1085  edge_val(n,2) = tri_x(n+1)
1086 
1087 end subroutine edge_values_implicit_h6
1088 
1089 end module regrid_edge_values
regrid_solvers
Solvers of linear systems.
Definition: regrid_solvers.F90:2
regrid_edge_values
Edge value estimation for high-order resconstruction.
Definition: regrid_edge_values.F90:2
polynomial_functions
Polynomial functions.
Definition: polynomial_functions.F90:2