Recon1d_EMPLM_CWK.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!> Piecewise Linear Method 1D reconstruction in index space and boundary extrapolation
6!!
7!! This implementation of PLM follows Colella and Woodward, 1984 \cite colella1984, except for assuming
8!! uniform resolution so that the method is independent of grid spacing. The cell-wise reconstructions
9!! are limited so that the edge values (which are also the extrema in a cell) are bounded by the neighbors.
10!! The slope of the first and last cells are set so that the first interior edge values match the interior
11!! cell (i.e. extrapolates from the interior).
13
14use recon1d_type, only : testing
15use recon1d_mplm_cwk, only : mplm_cwk
16
17implicit none ; private
18
19public emplm_cwk, testing
20
21!> PLM reconstruction following Colella and Woodward, 1984
22!!
23!! Implemented by extending recon1d_mplm_cwk.
24!!
25!! The source for the methods ultimately used by this class are:
26!! - init() -> recon1d_mplm_cwk -> recon1d_plm_cw.init()
27!! - reconstruct() *locally defined
28!! - average() -> recon1d_mplm_cwk -> recon1d_plm_cw.average()
29!! - f() -> recon1d_mplm_cwk -> recon1d_plm_cw.f()
30!! - dfdx() -> recon1d_mplm_cwk -> recon1d_plm_cw.dfdx()
31!! - check_reconstruction() -> recon1d_mplm_cwk.check_reconstruction()
32!! - unit_tests() *locally defined
33!! - destroy() -> recon1d_mplm_cwk -> recon1d_plm_cw.destroy()
34!! - remap_to_sub_grid() -> recon1d_type.remap_to_sub_grid()
35!! - init_parent() -> init()
36!! - reconstruct_parent() -> recon1d_mplm_cwk.reconstruct()
37type, extends (mplm_cwk) :: emplm_cwk
38
39contains
40 !> Implementation of the EMPLM_CWK reconstruction
41 procedure :: reconstruct => reconstruct
42 !> Implementation of unit tests for the EMPLM_CWK reconstruction
43 procedure :: unit_tests => unit_tests
44
45end type emplm_cwk
46
47contains
48
49!> Calculate a 1D PLM reconstructions based on h(:) and u(:)
50subroutine reconstruct(this, h, u)
51 class(emplm_cwk), intent(inout) :: this !< This reconstruction
52 real, intent(in) :: h(*) !< Grid spacing (thickness) [typically H]
53 real, intent(in) :: u(*) !< Cell mean values [A]
54 ! Local variables
55 real :: slp ! The PLM slopes (difference across cell) [A]
56 real :: sigma_l, sigma_c, sigma_r ! Left, central and right slope estimates as
57 ! differences across the cell [A]
58 real :: u_min, u_max ! Minimum and maximum value across cell [A]
59 real :: u_l, u_r, u_c ! Left, right, and center values [A]
60 real :: u_e(this%n+1) ! Average of edge values [A]
61 integer :: k, n
62
63 n = this%n
64
65 call this%reconstruct_parent(h, u)
66
67 this%ur(1) = this%ul(2)
68 this%ul(1) = u(1) + ( u(1) - this%ur(1) )
69
70 this%ul(n) = this%ur(n-1)
71 this%ur(n) = u(n) + ( u(n) - this%ul(n) )
72
73end subroutine reconstruct
74
75!> Runs PLM reconstruction unit tests and returns True for any fails, False otherwise
76logical function unit_tests(this, verbose, stdout, stderr)
77 class(emplm_cwk), intent(inout) :: this !< This reconstruction
78 logical, intent(in) :: verbose !< True, if verbose
79 integer, intent(in) :: stdout !< I/O channel for stdout
80 integer, intent(in) :: stderr !< I/O channel for stderr
81 ! Local variables
82 real, allocatable :: ul(:), ur(:), um(:) ! test values [A]
83 real, allocatable :: ull(:), urr(:) ! test values [A]
84 type(testing) :: test ! convenience functions
85 integer :: k
86
87 call test%set( stdout=stdout ) ! Sets the stdout channel in test
88 call test%set( stderr=stderr ) ! Sets the stderr channel in test
89 call test%set( verbose=verbose ) ! Sets the verbosity flag in test
90
91 call this%init(3)
92 call test%test( this%n /= 3, 'Setting number of levels')
93 allocate( um(3), ul(3), ur(3), ull(3), urr(3) )
94
95 call this%reconstruct( (/2.,2.,2./), (/1.,3.,5./) )
96 call test%real_arr(3, this%u_mean, (/1.,3.,5./), 'Setting cell values')
97
98 do k = 1, 3
99 ul(k) = this%f(k, 0.)
100 um(k) = this%f(k, 0.5)
101 ur(k) = this%f(k, 1.)
102 enddo
103 call test%real_arr(3, ul, (/0.,2.,4./), 'Evaluation on left edge')
104 call test%real_arr(3, um, (/1.,3.,5./), 'Evaluation in center')
105 call test%real_arr(3, ur, (/2.,4.,6./), 'Evaluation on right edge')
106
107 do k = 1, 3
108 ul(k) = this%dfdx(k, 0.)
109 um(k) = this%dfdx(k, 0.5)
110 ur(k) = this%dfdx(k, 1.)
111 enddo
112 call test%real_arr(3, ul, (/2.,2.,2./), 'dfdx on left edge')
113 call test%real_arr(3, um, (/2.,2.,2./), 'dfdx in center')
114 call test%real_arr(3, ur, (/2.,2.,2./), 'dfdx on right edge')
115
116 do k = 1, 3
117 um(k) = this%average(k, 0.25, 0.75) ! Average from x=0.25 to 0.75 in each cell
118 enddo
119 call test%real_arr(3, um, (/1.,3.,5./), 'Return interval average')
120
121 call this%destroy()
122 deallocate( um, ul, ur, ull, urr )
123
124 allocate( um(4), ul(4), ur(4) )
125 call this%init(4)
126
127 ! These values lead to non-monotonic reconstuctions which are
128 ! valid for transport problems but not always appropriate for
129 ! remapping to arbitrary resolution grids.
130 ! The O(h^2) slopes are -, 2, 2, - and the limited
131 ! slopes are 0, 1, 1, 0 so the everywhere the reconstructions
132 ! are bounded by neighbors but ur(2) and ul(3) are out-of-order.
133 call this%reconstruct( (/1.,1.,1.,1./), (/0.,3.,4.,7./) )
134 do k = 1, 4
135 ul(k) = this%f(k, 0.)
136 ur(k) = this%f(k, 1.)
137 enddo
138 call test%real_arr(4, ul, (/-2.5,2.5,3.5,4.5/), 'Evaluation on left edge')
139 call test%real_arr(4, ur, (/2.5,3.5,4.5,9.5/), 'Evaluation on right edge')
140
141 deallocate( um, ul, ur )
142
143 unit_tests = test%summarize('EMPLM_CWK:unit_tests')
144
145end function unit_tests
146
147!> \namespace recon1d_emplm_cwk
148!!
149
150end module recon1d_emplm_cwk