Neko 1.99.1
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
elementwise_filter.f90
Go to the documentation of this file.
1! Copyright (c) 2024, The Neko Authors
2! All rights reserved.
3!
4! Redistribution and use in source and binary forms, with or without
5! modification, are permitted provided that the following conditions
6! are met:
7!
8! * Redistributions of source code must retain the above copyright
9! notice, this list of conditions and the following disclaimer.
10!
11! * Redistributions in binary form must reproduce the above
12! copyright notice, this list of conditions and the following
13! disclaimer in the documentation and/or other materials provided
14! with the distribution.
15!
16! * Neither the name of the authors nor the names of its
17! contributors may be used to endorse or promote products derived
18! from this software without specific prior written permission.
19!
20! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
23! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
24! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
25! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
26! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
27! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
28! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
30! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31! POSSIBILITY OF SUCH DAMAGE.
32!
33!
36 use num_types, only : rp
37 use filter, only: filter_t
38 use math, only : rzero, rone, copy
39 use field, only : field_t
40 use coefs, only : coef_t
41 use utils, only : neko_error
43 use json_module, only : json_file
45 use speclib, only : zwgll, legendre_poly
46 use matrix, only : matrix_t
47 use mxm_wrapper, only : mxm
48 use tensor, only : tnsr3d, trsp
50 use device_math, only : device_cfill
51 use, intrinsic :: iso_c_binding, only : c_ptr, c_null_ptr, c_associated
52 implicit none
53 private
54
56 type, public, extends(filter_t) :: elementwise_filter_t
59 character(len=:), allocatable :: filter_type
61 integer :: nx
63 integer :: nt
65 real(kind=rp), allocatable :: fh(:,:), fht(:,:)
66 type(c_ptr) :: fh_d = c_null_ptr
67 type(c_ptr) :: fht_d = c_null_ptr
69 real(kind=rp), allocatable :: transfer(:)
70 contains
72 procedure, pass(this) :: init => elementwise_filter_init_from_json
74 procedure, pass(this) :: init_from_components => &
77 procedure, pass(this) :: free => elementwise_filter_free
79 procedure, pass(this) :: build_1d
81 procedure, pass(this) :: apply => elementwise_field_filter_3d
83
84contains
86 subroutine elementwise_filter_init_from_json(this, json, coef)
87 class(elementwise_filter_t), intent(inout) :: this
88 type(json_file), intent(inout) :: json
89 type(coef_t), intent(in) :: coef
90 real(kind=rp), allocatable :: transfer(:)
91 character(len=:), allocatable :: filter_type
92
93 ! Filter assumes lx = ly = lz
94 call this%init_base(json, coef)
95
96 call this%init_from_components(coef%dof%xh%lx)
97
98 call json_get_or_default(json, "filter.elementwise_filter_type", &
99 this%filter_type, "nonBoyd")
100
101 if (json%valid_path('filter.transfer_function')) then
102 call json_get(json, 'filter.transfer_function', transfer)
103 if (size(transfer) .eq. coef%dof%xh%lx) then
104 this%transfer = transfer
105 else
106 call neko_error("The transfer function of the elementwise " // &
107 "filter must correspond the order of the polynomial")
108 end if
109 call this%build_1d()
110 end if
111
113
117 class(elementwise_filter_t), intent(inout) :: this
118 integer :: nx
119
120 this%nx = nx
121 this%nt = nx ! initialize as if nothing is filtered yet
122
123 allocate(this%fh(nx, nx))
124 allocate(this%fht(nx, nx))
125 allocate(this%transfer(nx))
126
127 call rzero(this%fh, nx*nx)
128 call rzero(this%fht, nx*nx)
129 call rone(this%transfer, nx) ! initialize as if nothing is filtered yet
130
131 if (neko_bcknd_device .eq. 1) then
132 call device_map(this%fh, this%fh_d, this%nx * this%nx)
133 call device_map(this%fht, this%fht_d, this%nx * this%nx)
134 call device_cfill(this%fh_d, 0.0_rp, this%nx * this%nx)
135 call device_cfill(this%fht_d, 0.0_rp, this%nx * this%nx)
136 end if
137
139
141 subroutine elementwise_filter_free(this)
142 class(elementwise_filter_t), intent(inout) :: this
143
144 if (allocated(this%fh)) then
145 deallocate(this%fh)
146 end if
147
148 if (allocated(this%fht)) then
149 deallocate(this%fht)
150 end if
151
152 if (allocated(this%transfer)) then
153 deallocate(this%transfer)
154 end if
155
156 if (c_associated(this%fh_d)) then
157 call device_free(this%fh_d)
158 end if
159
160 if (c_associated(this%fht_d)) then
161 call device_free(this%fht_d)
162 end if
163
164 this%filter_type = ""
165 this%nx = 0
166 this%nt = 0
167
168 call this%free_base()
169
170 end subroutine elementwise_filter_free
171
173 subroutine build_1d(this)
174 class(elementwise_filter_t), intent(inout) :: this
175
176 call build_1d_cpu(this%fh, this%fht, this%transfer, &
177 this%nx, this%filter_type)
178 if (neko_bcknd_device .eq. 1) then
179 call device_memcpy(this%fh, this%fh_d, &
180 this%nx * this%nx, host_to_device, sync = .false.)
181 call device_memcpy(this%fht, this%fht_d, &
182 this%nx * this%nx, host_to_device, sync = .false.)
183 end if
184
185 end subroutine build_1d
186
188 subroutine elementwise_field_filter_3d(this, F_out, F_in)
189 class(elementwise_filter_t), intent(inout) :: this
190 type(field_t), intent(inout) :: F_out
191 type(field_t), intent(in) :: F_in
192
193 ! F_out = fh x fh x fh x F_in
194 call tnsr3d(f_out%x, this%nx, f_in%x, this%nx, this%fh, this%fht, this%fht, &
195 this%coef%msh%nelv)
196
197 end subroutine elementwise_field_filter_3d
198
206 subroutine build_1d_cpu(fh, fht, transfer, nx, filter_type)
207 integer, intent(in) :: nx
208 real(kind=rp), intent(inout) :: fh(nx, nx), fht(nx, nx)
209 real(kind=rp), intent(in) :: transfer(nx)
210 real(kind=rp) :: diag(nx, nx), rmult(nx), lj(nx), zpts(nx)
211 type(matrix_t) :: phi, pht
212 integer :: n, i, j, k
213 real(kind=rp) :: z
214 character(len=*), intent(in) :: filter_type
215
216 call phi%init(nx, nx)
217 call pht%init(nx, nx)
218
219 call zwgll(zpts, rmult, nx)
220
221 n = nx-1
222 do j = 1, nx
223 z = zpts(j)
224 call legendre_poly(lj, z, n)
225 select case (filter_type)
226 case("Boyd")
227 pht%x(1,j) = lj(1)
228 pht%x(2,j) = lj(2)
229 do k=3,nx
230 pht%x(k,j) = lj(k)-lj(k-2)
231 end do
232 case("nonBoyd")
233 pht%x(:,j) = lj
234 end select
235 end do
236
237 call trsp(phi%x, nx, pht%x, nx)
238 pht%x = phi%x
239
240 call pht%inverse(0) ! "0" for cpu implementation
241
242 diag = 0.0_rp
243
244 do i=1,nx
245 diag(i,i) = transfer(i)
246 end do
247
248 call mxm (diag, nx, pht%x, nx, fh, nx) ! -1
249 call mxm (phi%x, nx, fh, nx, pht%x, nx) ! V D V
250
251 call copy (fh, pht%x, nx*nx)
252 call trsp (fht, nx, fh, nx)
253
254 end subroutine build_1d_cpu
255
256end module elementwise_filter
Map a Fortran array to a device (allocate and associate)
Definition device.F90:72
Copy data between host and device (or device and device)
Definition device.F90:66
Retrieves a parameter by name or assigns a provided default value. In the latter case also adds the m...
Retrieves a parameter by name or throws an error.
Coefficients.
Definition coef.f90:34
subroutine, public device_cfill(a_d, c, n, strm)
Set all elements to a constant c .
Device abstraction, common interface for various accelerators.
Definition device.F90:34
integer, parameter, public host_to_device
Definition device.F90:47
subroutine, public device_free(x_d)
Deallocate memory on the device.
Definition device.F90:214
Implements elementwise_filter_t.
subroutine build_1d_cpu(fh, fht, transfer, nx, filter_type)
Build the 1d filter for an element on the CPU. Suppose field x is filtered into x_hat by x_hat = fh*x...
subroutine elementwise_filter_init_from_json(this, json, coef)
Constructor.
subroutine elementwise_filter_free(this)
Destructor.
subroutine build_1d(this)
Build the 1d filter for an element.
subroutine elementwise_filter_init_from_components(this, nx)
Actual Constructor.
subroutine elementwise_field_filter_3d(this, f_out, f_in)
Filter a 3D field.
Defines a field.
Definition field.f90:34
Filter to be applied to a scalar field.
Definition filter.f90:38
Utilities for retrieving parameters from the case files.
Definition math.f90:60
subroutine, public rone(a, n)
Set all elements to one.
Definition math.f90:244
subroutine, public copy(a, b, n)
Copy a vector .
Definition math.f90:255
subroutine, public rzero(a, n)
Zero a real vector.
Definition math.f90:211
Defines a matrix.
Definition matrix.f90:34
Wrapper for all matrix-matrix product implementations.
subroutine, public mxm(a, n1, b, n2, c, n3)
Compute matrix-matrix product for contiguously packed matrices A,B, and C.
Build configurations.
integer, parameter neko_bcknd_device
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12
LIBRARY ROUTINES FOR SPECTRAL METHODS.
Definition speclib.f90:148
subroutine zwgll(z, w, np)
Generate NP Gauss-Lobatto Legendre points (Z) and weights (W) associated with Jacobi polynomial P(N)(...
Definition speclib.f90:179
subroutine legendre_poly(l, x, n)
Evaluate Legendre polynomials of degrees 0-N at point x and store in array L.
Definition speclib.f90:1000
Tensor operations.
Definition tensor.f90:61
subroutine, public trsp(a, lda, b, ldb)
Transpose of a rectangular tensor .
Definition tensor.f90:121
subroutine, public tnsr3d(v, nv, u, nu, a, bt, ct, nelv)
Tensor product performed on nelv elements.
Definition tensor.f90:225
Utilities.
Definition utils.f90:35
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:55
Implements the elementwise filter for SEM.
Base abstract class for filter.
Definition filter.f90:48