Neko 1.99.1
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
cpr.f90
Go to the documentation of this file.
1! Copyright (c) 2022, 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!
34module cpr
35 use num_types, only : rp
36 use field, only : field_t
37 use space, only : space_t
38 use math, only : copy, sub3, flipv, swap, reord, sort
39 use tensor, only : trsp1
40 use mesh, only : mesh_t
41 use coefs, only : coef_t
42 use logger, only : neko_log, log_size
43 use dofmap, only : dofmap_t
44 use mxm_wrapper, only : mxm
45 implicit none
46 private
47
49 type, public :: cpr_t
50 real(kind=rp), allocatable :: v(:,:)
51
52 real(kind=rp), allocatable :: vt(:,:)
53 real(kind=rp), allocatable :: vinv(:,:)
54 real(kind=rp), allocatable :: vinvt(:,:)
56 real(kind=rp), allocatable :: w(:,:)
57
58 real(kind=rp), allocatable :: fldhat(:,:,:,:)
59
60 type(field_t), pointer :: fld => null()
61 type(space_t), pointer :: xh => null()
62 type(mesh_t), pointer :: msh => null()
63 type(dofmap_t), pointer :: dof => null()
64
65 end type cpr_t
66
67 interface cpr_init
68 module procedure cpr_init_all
69 end interface cpr_init
70
71 public :: cpr_init, cpr_free
72
73contains
74
76 subroutine cpr_init_all(cpr, u)
77 type(cpr_t), intent(inout) :: cpr
78 type(field_t), intent(in), target :: u
79
80 call cpr_free(cpr)
81
82 cpr%fld => u
83 cpr%msh => u%msh
84 cpr%Xh => u%Xh
85 cpr%dof => u%dof
86
87
88 ! Allocate arrays
89 allocate(cpr%v(cpr%Xh%lx, cpr%Xh%lx))
90 allocate(cpr%vt(cpr%Xh%lx, cpr%Xh%lx))
91 allocate(cpr%vinv(cpr%Xh%lx, cpr%Xh%lx))
92 allocate(cpr%vinvt(cpr%Xh%lx, cpr%Xh%lx))
93 allocate(cpr%w(cpr%Xh%lx, cpr%Xh%lx))
94 allocate(cpr%fldhat(cpr%Xh%lx, cpr%Xh%ly, cpr%Xh%lz, cpr%msh%nelv))
95
96 ! Initialize all the matrices
98
99 ! Generate the uhat field (legendre coeff)
100 call cpr_goto_space(cpr,'spec')
101
102 end subroutine cpr_init_all
103
104
106 subroutine cpr_free(cpr)
107 type(cpr_t), intent(inout) :: cpr
108
109 if(allocated(cpr%v)) then
110 deallocate(cpr%v)
111 end if
112
113 if(allocated(cpr%vt)) then
114 deallocate(cpr%vt)
115 end if
116
117 if(allocated(cpr%vinv)) then
118 deallocate(cpr%vinv)
119 end if
120
121 if(allocated(cpr%vinvt)) then
122 deallocate(cpr%vinvt)
123 end if
124
125 if(allocated(cpr%w)) then
126 deallocate(cpr%w)
127 end if
128
129 if(allocated(cpr%fldhat)) then
130 deallocate(cpr%fldhat)
131 end if
132
133 nullify(cpr%fld)
134 nullify(cpr%msh)
135 nullify(cpr%Xh)
136 nullify(cpr%dof)
137
138
139 end subroutine cpr_free
140
141
143 subroutine cpr_generate_specmat(cpr)
144 type(cpr_t), intent(inout) :: cpr
145 real(kind=rp) :: l(0:cpr%Xh%lx-1)
146 real(kind=rp) :: delta(cpr%Xh%lx)
147 integer :: i, kj, j, j2, kk
148
149 associate(xh => cpr%Xh, v=> cpr%v, vt => cpr%vt, &
150 vinv => cpr%vinv, vinvt => cpr%vinvt, w => cpr%w)
151 ! Get the Legendre polynomials for each point
152 ! Then proceed to compose the transform matrix
153 kj = 0
154 do j = 1, xh%lx
155 l(0) = 1.
156 l(1) = xh%zg(j,1)
157 do j2 = 2, xh%lx-1
158 l(j2) = ( (2*j2-1) * xh%zg(j,1) * l(j2-1) &
159 - (j2-1) * l(j2-2) ) / j2
160 end do
161 do kk = 1, xh%lx
162 kj = kj+1
163 v(kj,1) = l(kk-1)
164 end do
165 end do
166
167 ! transpose the matrix
168 call trsp1(v, xh%lx)
169
170 ! Calculate the nominal scaling factors
171 do i = 1, xh%lx
172 delta(i) = 2.0_rp / (2*(i-1)+1)
173 end do
174 ! modify last entry
175 delta(xh%lx) = 2.0_rp / (xh%lx-1)
176
177 ! calculate the inverse to multiply the matrix
178 do i = 1, xh%lx
179 delta(i) = sqrt(1.0_rp / delta(i))
180 end do
181 ! scale the matrix
182 do i = 1, xh%lx
183 do j = 1, xh%lx
184 v(i,j) = v(i,j) * delta(j) ! orthogonal wrt weights
185 end do
186 end do
187
188 ! get the trasposed
189 call copy(vt, v, xh%lx * xh%lx)
190 call trsp1(vt, xh%lx)
191
192 !populate the mass matrix
193 kk = 1
194 do i = 1, xh%lx
195 do j = 1, xh%lx
196 if (i .eq. j) then
197 w(i,j) = xh%wx(kk)
198 kk = kk+1
199 else
200 cpr%w(i,j) = 0
201 end if
202 end do
203 end do
204
205 !Get the inverse of the transform matrix
206 call mxm(vt, xh%lx, w, xh%lx, vinv, xh%lx)
207
208 !get the transposed of the inverse
209 call copy(vinvt, vinv, xh%lx * xh%lx)
210 call trsp1(vinvt, xh%lx)
211 end associate
212
213 end subroutine cpr_generate_specmat
214
215
217 !the result of the transform is given in fldhat
218 subroutine cpr_goto_space(cpr, space)
219 type(cpr_t), intent(inout) :: cpr
220 real(kind=rp) :: w2(cpr%Xh%lx,cpr%Xh%lx,cpr%Xh%lx)
221 real(kind=rp) :: w1(cpr%Xh%lx,cpr%Xh%lx,cpr%Xh%lx)
222 real(kind=rp) :: specmat(cpr%Xh%lx,cpr%Xh%lx)
223 real(kind=rp) :: specmatt(cpr%Xh%lx,cpr%Xh%lx)
224 integer :: k, e, nxyz, nelv
225 character(len=4) :: space
226
227 ! define some constants
228 nxyz = cpr%Xh%lx*cpr%Xh%lx*cpr%Xh%lx
229 nelv = cpr%msh%nelv
230
231 ! Define the matrix according to which transform to do
232 if (space .eq. 'spec') then
233 call copy(specmat, cpr%vinv, cpr%Xh%lx*cpr%Xh%lx)
234 call copy(specmatt, cpr%vinvt, cpr%Xh%lx*cpr%Xh%lx)
235 endif
236 if (space .eq. 'phys') then
237 call copy(specmat, cpr%v, cpr%Xh%lx*cpr%Xh%lx)
238 call copy(specmatt, cpr%vt, cpr%Xh%lx*cpr%Xh%lx)
239 endif
240
241 ! Apply the operator (transform to given space)
242 do e=1,nelv
243 if (space .eq. 'spec') then
244 call copy(w2, cpr%fld%x(1,1,1,e), nxyz) ! start from phys field
245 else
246 call copy(w2, cpr%fldhat(1,1,1,e), nxyz) ! start from spec coeff
247 endif
248 ! apply the operator to the x direction, result in w1
249 call mxm(specmat, cpr%Xh%lx, w2, cpr%Xh%lx, w1, cpr%Xh%lx*cpr%Xh%lx)
250 ! apply matrix to y direction, result in w2
251 do k=1,cpr%Xh%lx
252 call mxm(w1(1,1,k),cpr%Xh%lx,specmatt,cpr%Xh%lx,w2(1,1,k),cpr%Xh%lx)
253 enddo
254 ! apply matrix to z direction, result always in fldhat
255 call mxm (w2, cpr%Xh%lx*cpr%Xh%lx, specmatt, cpr%Xh%lx,&
256 cpr%fldhat(1,1,1,e), cpr%Xh%lx)
257 enddo
258
259 end subroutine cpr_goto_space
260
262 subroutine cpr_truncate_wn(cpr, coef)
263 type(cpr_t), intent(inout) :: cpr
264 type(coef_t), intent(inout) :: coef
265 real(kind=rp) :: w2(cpr%Xh%lx, cpr%Xh%lx, cpr%Xh%lx)
266 real(kind=rp) :: w1(cpr%Xh%lx, cpr%Xh%lx, cpr%Xh%lx)
267 real(kind=rp) :: vsort(cpr%Xh%lx * cpr%Xh%lx * cpr%Xh%lx)
268 real(kind=rp) :: vtrunc(cpr%Xh%lx, cpr%Xh%lx, cpr%Xh%lx)
269 real(kind=rp) :: vtemp(cpr%Xh%lx * cpr%Xh%lx * cpr%Xh%lx)
270 real(kind=rp) :: errvec(cpr%Xh%lx, cpr%Xh%lx, cpr%Xh%lx)
271 real(kind=rp) :: fx(cpr%Xh%lx, cpr%Xh%lx)
272 real(kind=rp) :: fy(cpr%Xh%lx, cpr%Xh%lx)
273 real(kind=rp) :: fz(cpr%Xh%lx, cpr%Xh%lx)
274 real(kind=rp) :: l2norm, oldl2norm, targeterr
275 integer :: isort(cpr%Xh%lx * cpr%Xh%lx * cpr%Xh%lx)
276 integer :: i, j, k, e, nxyz, nelv
277 integer :: kut, kutx, kuty, kutz, nx
278 character(len=LOG_SIZE) :: log_buf
279
280 ! define some constants
281 nx = cpr%Xh%lx
282 nxyz = cpr%Xh%lx*cpr%Xh%lx*cpr%Xh%lx
283 nelv = cpr%msh%nelv
284 targeterr = 1e-3_rp
285
286 !truncate for every element
287 do e = 1, nelv
288 ! create temp vector where the trunc coeff will be
289 call copy(vtemp, cpr%fldhat(1,1,1,e), nxyz)
290 call copy(vtrunc, cpr%fldhat(1,1,1,e), nxyz)
291 ! sort the coefficients by absolute value
292 call sortcoeff(vsort, cpr%fldhat(1,1,1,e), isort, nxyz)
293 ! initialize values for iterative procedure
294 l2norm = 0.0_rp
295 kut = 0
296
297 do while (l2norm .le. targeterr .and. kut .le. (cpr%Xh%lx-1)*3)
298 ! save value from prev it since it met requirements to enter
299 call copy(vtrunc, vtemp, nxyz)
300 oldl2norm = l2norm
301
302 ! advance kut
303 kut = kut+1
304
305 ! create filters
306 call build_filter_tf(fx, fy, fz, kut, cpr%Xh%lx)
307
308 ! truncate, remember that transposed of diag mat is the same
309 call copy(w2, vsort, nxyz)
310 ! apply the operator to the x direction, result in w1
311 call mxm(fx, nx, w2, nx, w1, nx*nx)
312 ! apply matrix to y direction, result in w2
313 do k=1,nx
314 call mxm(w1(1,1,k), nx, fy, nx, w2(1,1,k), nx)
315 enddo
316 ! apply matrix to z direction, result in vtemp
317 call mxm (w2, nx*nx, fz, nx, vtemp, nx)
318
319 ! vtemp is sorted, so bring it back to real order
320 ! by inverting the swap operation in sortcoeff
321 call reord(vtemp, isort, nxyz)
322
323 ! calculate the error vector
324 call sub3(errvec, cpr%fldhat(1,1,1,e), vtemp, nxyz)
325
326 ! get the norm of the error
327 l2norm = get_elem_l2norm(errvec, coef, 'spec' ,e)
328
329 end do
330
331 ! copy the truncated field back to the main object
332 call copy(cpr%fldhat(1,1,1,e), vtrunc, nxyz)
333
334
335 !================debugging info
336
337 !just to check that all is good, resort vsort
338 if (e .eq. 50) then
339 write(log_buf, '(A)') &
340 'Debugging info for e=50. Do not forget to delete'
341 call neko_log%message(log_buf)
342
343 call reord(vsort,isort,nxyz)
344
345 write(log_buf, '(A)') &
346 'Norm'
347 call neko_log%message(log_buf)
348 !chech that the copy is fine in one entry
349 write(log_buf, '(A,E15.7)') &
350 'l2norm:', oldl2norm
351 call neko_log%message(log_buf)
352 write(log_buf, '(A)') &
353 'cut off level'
354 call neko_log%message(log_buf)
355 !chech that the copy is fine in one entry
356 write(log_buf, '(A,I5)') &
357 'kut:', kut
358 call neko_log%message(log_buf)
359 write(log_buf, '(A)') &
360 'spectral coefficients'
361 call neko_log%message(log_buf)
362 do i = 1, 10
363 write(log_buf, '(A,E15.7,A,E15.7)') &
364 'org coeff:', vsort(i), ' truncated coeff', &
365 cpr%fldhat(i,1,1,e)
366 call neko_log%message(log_buf)
367 end do
368 end if
369 end do
370
371 end subroutine cpr_truncate_wn
372
375 subroutine sortcoeff(vsort, v, isort, nxyz)
376 integer, intent(in) :: nxyz
377 real(kind=rp), intent(inout) :: vsort(nxyz)
378 real(kind=rp), intent(inout) :: v(nxyz)
379 integer, intent(inout) :: isort(nxyz)
380 integer :: i
381
382 ! copy absolute values to sort by magnitude
383 do i =1, nxyz
384 vsort(i) = abs(v(i))
385 isort(i) = i
386 end do
387
388 ! sort the absolute values of the vectors, here the index is the
389 ! important part (isort)
390 call sort(vsort, isort, nxyz)
391
392 ! Flip the indices so they are in a descending order
393 call flipv(vsort, isort, nxyz)
394
395 ! now re order the fields (non abs value) based on the indices
396 call copy(vsort, v, nxyz)
397 call swap(vsort, isort, nxyz)
398
399 end subroutine sortcoeff
400
402 subroutine build_filter_tf(fx, fy, fz, kut, lx)
403 integer, intent(in) :: lx
404 integer, intent(in) :: kut
405 real(kind=rp), intent(inout) :: fx(lx,lx)
406 real(kind=rp), intent(inout) :: fy(lx,lx)
407 real(kind=rp), intent(inout) :: fz(lx,lx)
408 integer :: i, j, k, k0, kk, kutx, kuty, kutz
409
410 ! set the values acording to kut
411 if (kut .eq. 0) then
412 kutz = 0
413 kuty = 0
414 kutx = 0
415 end if
416 if (kut .gt. 0 .and. kut .le. (lx-1)) then
417 kutz = kut
418 kuty = 0
419 kutx = 0
420 end if
421 if (kut .gt. (lx-1) .and. kut .le. (lx-1)*2) then
422 kutz = lx-1
423 kuty = kut-(lx-1)
424 kutx = 0
425 end if
426 if (kut .gt. (lx-1)*2 .and. kut .le. (lx-1)*3) then
427 kutz = lx-1
428 kuty = lx-1
429 kutx = kut-(lx-1)*2
430 end if
431
432 ! create identity matrices
433 do i = 1, lx
434 do j = 1, lx
435 if (i .eq. j) then
436 fx(i,j) = 1
437 fy(i,j) = 1
438 fz(i,j) = 1
439 else
440 fx(i,j) = 0
441 fy(i,j) = 0
442 fz(i,j) = 0
443 end if
444 end do
445 end do
446
447 ! truncate the matrices acording to input kut
448 k0 = lx-kutx
449 do k = k0+1, lx
450 kk = k+lx*(k-1)
451 fx(kk,1) = 0
452 end do
453 k0 = lx-kuty
454 do k = k0+1, lx
455 kk = k+lx*(k-1)
456 fy(kk,1) = 0
457 end do
458 k0 = lx-kutz
459 do k = k0+1, lx
460 kk = k+lx*(k-1)
461 fz(kk,1) = 0
462 end do
463
464 end subroutine build_filter_tf
465
466
467 function get_elem_l2norm(elemdata, coef, space, e) result(l2e)
468 type(coef_t), intent(in) :: coef
469 real(kind=rp) :: elemdata(coef%Xh%lx, coef%Xh%lx, coef%Xh%lx)
470 real(kind=rp) :: vole, suma, l2e
471 integer i, e, nxyz
472 character(len=4) :: space
473
474 ! Get the volume of the element
475 nxyz = coef%Xh%lx*coef%Xh%lx*coef%Xh%lx
476 vole=0
477 do i = 1, nxyz
478 vole = vole+coef%B(i,1,1,e)
479 end do
480
481 ! Get the weighted l2 norm of the element
482 suma = 0
483 if (space .eq. 'spec') then
484 do i = 1, nxyz
485 suma = suma+elemdata(i,1,1)*elemdata(i,1,1)*coef%jac(i,1,1,e)
486 end do
487 end if
488 if (space .eq. 'phys') then
489 do i = 1, nxyz
490 suma = suma+elemdata(i,1,1)*elemdata(i,1,1)*coef%B(i,1,1,e)
491 end do
492 end if
493
494 l2e = sqrt(suma)/sqrt(vole)
495
496 end function get_elem_l2norm
497
498
499end module cpr
Coefficients.
Definition coef.f90:34
Compression.
Definition cpr.f90:34
subroutine build_filter_tf(fx, fy, fz, kut, lx)
create filter transfer function
Definition cpr.f90:403
real(kind=rp) function get_elem_l2norm(elemdata, coef, space, e)
Definition cpr.f90:468
subroutine cpr_truncate_wn(cpr, coef)
Truncate the frequencies.
Definition cpr.f90:263
subroutine cpr_init_all(cpr, u)
Initialize cpr.
Definition cpr.f90:77
subroutine sortcoeff(vsort, v, isort, nxyz)
Sort the spectral coefficient in descending order array vsort. The original indices are stored in the...
Definition cpr.f90:376
subroutine cpr_goto_space(cpr, space)
Tranform to spectral space (using tensor product)
Definition cpr.f90:219
subroutine cpr_generate_specmat(cpr)
Generate spectral tranform matrices.
Definition cpr.f90:144
subroutine, public cpr_free(cpr)
Deallocate coefficients.
Definition cpr.f90:107
Defines a mapping of the degrees of freedom.
Definition dofmap.f90:35
Defines a field.
Definition field.f90:34
Logging routines.
Definition log.f90:34
type(log_t), public neko_log
Global log stream.
Definition log.f90:70
integer, parameter, public log_size
Definition log.f90:46
Definition math.f90:60
subroutine, public sub3(a, b, c, n)
Vector subtraction .
Definition math.f90:787
subroutine, public copy(a, b, n)
Copy a vector .
Definition math.f90:255
Defines a mesh.
Definition mesh.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.
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12
Defines a function space.
Definition space.f90:34
Tensor operations.
Definition tensor.f90:61
subroutine, public trsp1(a, n)
In-place transpose of a square tensor.
Definition tensor.f90:137
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:55
include information needed for compressing fields
Definition cpr.f90:49
The function space for the SEM solution fields.
Definition space.f90:62