Recon1d_PCM.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!> 1D reconstructions using the Piecewise Constant Method (PCM)
6module recon1d_pcm
7
8use recon1d_type, only : recon1d, testing
9
10implicit none ; private
11
12public pcm
13
14!> PCM (piecewise constant) reconstruction
15!!
16!! The source for the methods ultimately used by this class are:
17!! init() *locally defined
18!! reconstruct() *locally defined
19!! average() *locally defined
20!! f() *locally defined
21!! dfdx() *locally defined
22!! - x() *locally defined
23!! check_reconstruction() *locally defined
24!! unit_tests() *locally defined
25!! destroy() *locally defined
26!! remap_to_sub_grid() -> Recon1d%remap_to_sub_grid()
27!! init_parent() -> init()
28!! reconstruct_parent() -> parent()
30
31contains
32 !> Implementation of the PCM initialization
34 !> Implementation of the PCM reconstruction
35 procedure :: reconstruct => reconstruct
36 !> Implementation of the PCM average over an interval [A]
38 !> Implementation of evaluating the PCM reconstruction at a point [A]
40 !> Implementation of the derivative of the PCM reconstruction at a point [A]
42 !> Implementation of solver for x: f(x)=t
44 !> Implementation of deallocation for PCM
46 !> Implementation of check reconstruction for the PCM reconstruction
47 procedure :: check_reconstruction => check_reconstruction
48 !> Implementation of unit tests for the PCM reconstruction
49 procedure :: unit_tests => unit_tests
50
51 !> Duplicate interface to init()
52 procedure :: init_parent => init
53 !> Duplicate interface to reconstruct()
54 procedure :: reconstruct_parent => reconstruct
55
56end type pcm
57
58contains
59
60!> Initialize a 1D PCM reconstruction for n cells
61subroutine init(this, n, h_neglect, check)
62 class(pcm), intent(out) :: this !< This reconstruction
63 integer, intent(in) :: n !< Number of cells in this column
64 real, optional, intent(in) :: h_neglect !< A negligibly small width used in cell reconstructions [H].
65 !! Not used by PCM.
66 logical, optional, intent(in) :: check !< If true, enable some consistency checking
67
68 if (present(h_neglect)) this%n = n ! no-op to avoid compiler warning about unused dummy argument
69 if (present(check)) this%check = check
70
71 this%n = n
72
73 allocate( this%u_mean(n) )
74
75end subroutine init
76
77!> Calculate a 1D PCM reconstructions based on h(:) and u(:)
78subroutine reconstruct(this, h, u)
79 class(pcm), intent(inout) :: this !< This reconstruction
80 real, intent(in) :: h(*) !< Grid spacing (thickness) [typically H]
81 real, intent(in) :: u(*) !< Cell mean values [A]
82 ! Local variables
83 integer :: k
84
85 this%u_mean(1) = h(1) ! no-op to avoid compiler warning about unused dummy argument
86
87 do k = 1, this%n
88 this%u_mean(k) = u(k)
89 enddo
90
91end subroutine reconstruct
92
93!> Value of PCM reconstruction at a point in cell k [A]
94real function f(this, k, x)
95 class(pcm), intent(in) :: this !< This reconstruction
96 integer, intent(in) :: k !< Cell number
97 real, intent(in) :: x !< Non-dimensional position within element [nondim]
98
99 f = this%u_mean(k)
100
101end function f
102
103!> Derivative of PCM reconstruction at a point in cell k [A]
104real function dfdx(this, k, x)
105 class(pcm), intent(in) :: this !< This reconstruction
106 integer, intent(in) :: k !< Cell number
107 real, intent(in) :: x !< Non-dimensional position within element [nondim]
108
109 dfdx = 0.
110
111end function dfdx
112
113!> Solver for x: f(x)=t
114real function x(this, k, t)
115 class(pcm), intent(in) :: this !< This reconstruction
116 integer, intent(in) :: k !< Cell number
117 real, intent(in) :: t !< Value to solve for [A]
118 real :: slp ! Difference across cell [A]
119
120 slp = this%u_mean(min(k+1,this%n)) - this%u_mean(max(k-1,1))
121 if ( abs(slp) > 0. ) slp = sign(1., slp)
122 x = 0.5 ! Fall back if t==u_mean
123 ! if t>u_mean & slp=1 then x=1
124 ! if t<u_mean & slp=1 then x=0
125 ! if t>u_mean & slp=-1 then x=0
126 ! if t<u_mean & slp=-1 then x=1
127 ! if slp=0 then x=0.5
128 if ( abs(t - this%u_mean(k)) > 0. ) x = 0.5 + slp * sign(0.5, t - this%u_mean(k))
129end function x
130
131!> Average between xa and xb for cell k of a 1D PCM reconstruction [A]
132real function average(this, k, xa, xb)
133 class(pcm), intent(in) :: this !< This reconstruction
134 integer, intent(in) :: k !< Cell number
135 real, intent(in) :: xa !< Start of averaging interval on element (0 to 1)
136 real, intent(in) :: xb !< End of averaging interval on element (0 to 1)
137
138 average = xb + xa ! no-op to avoid compiler warnings about unused dummy argument
139 average = this%u_mean(k)
140
141end function average
142
143!> Deallocate the PCM reconstruction
144subroutine destroy(this)
145 class(pcm), intent(inout) :: this !< This reconstruction
146
147 deallocate( this%u_mean )
148
149end subroutine destroy
150
151!> Checks the PCM reconstruction for consistency
152logical function check_reconstruction(this, h, u)
153 class(pcm), intent(in) :: this !< This reconstruction
154 real, intent(in) :: h(*) !< Grid spacing (thickness) [typically H]
155 real, intent(in) :: u(*) !< Cell mean values [A]
156 ! Local variables
157 integer :: k
158
159 check_reconstruction = .false.
160
161 do k = 1, this%n
162 if ( abs( this%u_mean(k) - u(k) ) > 0. ) check_reconstruction = .true.
163 enddo
164
165end function check_reconstruction
166
167!> Runs PCM reconstruction unit tests and returns True for any fails, False otherwise
168logical function unit_tests(this, verbose, stdout, stderr)
169 class(pcm), intent(inout) :: this !< This reconstruction
170 logical, intent(in) :: verbose !< True, if verbose
171 integer, intent(in) :: stdout !< I/O channel for stdout
172 integer, intent(in) :: stderr !< I/O channel for stderr
173 ! Local variables
174 real, allocatable :: ul(:), ur(:), um(:) ! test values [A]
175 type(testing) :: test ! convenience functions
176 integer :: k
177
178 call test%set( stdout=stdout ) ! Sets the stdout channel in test
179 call test%set( stderr=stderr ) ! Sets the stderr channel in test
180 call test%set( verbose=verbose ) ! Sets the verbosity flag in test
181
182 call this%init(3)
183 call test%test( this%n /= 3, 'Setting number of levels')
184 allocate( um(3), ul(3), ur(3) )
185
186 call this%reconstruct( (/2.,2.,2./), (/1.,3.,5./) )
187 call test%real_arr(3, this%u_mean, (/1.,3.,5./), 'Setting cell values')
188
189 do k = 1, 3
190 ul(k) = this%f(k, 0.)
191 um(k) = this%f(k, 0.5)
192 ur(k) = this%f(k, 1.)
193 enddo
194 call test%real_arr(3, ul, (/1.,3.,5./), 'Evaluation on left edge')
195 call test%real_arr(3, um, (/1.,3.,5./), 'Evaluation in center')
196 call test%real_arr(3, ur, (/1.,3.,5./), 'Evaluation on right edge')
197
198 do k = 1, 3
199 ul(k) = this%dfdx(k, 0.)
200 um(k) = this%dfdx(k, 0.5)
201 ur(k) = this%dfdx(k, 1.)
202 enddo
203 call test%real_arr(3, ul, (/0.,0.,0./), 'dfdx on left edge')
204 call test%real_arr(3, um, (/0.,0.,0./), 'dfdx in center')
205 call test%real_arr(3, ur, (/0.,0.,0./), 'dfdx on right edge')
206
207 call test%real_scalar( this%x(1,0.), 0., 'f-1(1,0)=0')
208 call test%real_scalar( this%x(1,1.), 0.5, 'f-1(1,1)=0.5')
209 call test%real_scalar( this%x(1,3.), 1., 'f-1(1,3)=1')
210 call test%real_scalar( this%x(2,1.), 0., 'f-1(2,1)=0')
211 call test%real_scalar( this%x(2,3.), 0.5, 'f-1(2,3)=0.5')
212 call test%real_scalar( this%x(2,5.), 1., 'f-1(2,5)=1')
213 call test%real_scalar( this%x(3,3.), 0., 'f-1(3,3)=0')
214 call test%real_scalar( this%x(3,5.), 0.5, 'f-1(3,5)=0.5')
215 call test%real_scalar( this%x(3,7.), 1., 'f-1(3,7)=1')
216
217 do k = 1, 3
218 um(k) = this%average(k, 0.5, 0.75)
219 enddo
220 call test%real_arr(3, um, (/1.,3.,5./), 'Return interval average')
221
222 unit_tests = test%summarize('PCM:unit_tests')
223
224end function unit_tests
225
226!> \namespace recon1d_pcm
227!!
228
229end module recon1d_pcm