52 real(kind=
rp),
allocatable :: v(:,:)
54 real(kind=
rp),
allocatable :: vt(:,:)
55 real(kind=
rp),
allocatable :: vinv(:,:)
56 real(kind=
rp),
allocatable :: vinvt(:,:)
58 real(kind=
rp),
allocatable :: w(:,:)
60 real(kind=
rp),
allocatable :: fldhat(:,:,:,:)
64 type(
mesh_t),
pointer :: msh => null()
80 type(
field_t),
intent(in),
target :: u
111 if(
allocated(
cpr%v))
then
115 if(
allocated(
cpr%vt))
then
119 if(
allocated(
cpr%vinv))
then
123 if(
allocated(
cpr%vinvt))
then
124 deallocate(
cpr%vinvt)
127 if(
allocated(
cpr%w))
then
131 if(
allocated(
cpr%fldhat))
then
132 deallocate(
cpr%fldhat)
146 type(
cpr_t),
intent(inout) :: cpr
147 real(kind=
rp) :: l(0:
cpr%Xh%lx-1)
148 real(kind=
rp) :: delta(
cpr%Xh%lx)
149 integer :: i, kj, j, j2, kk
150 character(len=LOG_SIZE) :: log_buf
152 associate(xh =>
cpr%Xh, v=>
cpr%v, vt =>
cpr%vt, &
153 vinv =>
cpr%vinv, vinvt =>
cpr%vinvt, w =>
cpr%w)
161 l(j2) = ( (2*j2-1) * xh%zg(j,1) * l(j2-1) &
162 - (j2-1) * l(j2-2) ) / j2
175 delta(i) = 2.0_rp / (2*(i-1)+1)
178 delta(xh%lx) = 2.0_rp / (xh%lx-1)
182 delta(i) = sqrt(1.0_rp / delta(i))
187 v(i,j) = v(i,j) * delta(j)
192 call copy(vt, v, xh%lx * xh%lx)
193 call trsp1(vt, xh%lx)
209 call mxm(vt, xh%lx, w, xh%lx, vinv, xh%lx)
212 call copy(vinvt, vinv, xh%lx * xh%lx)
213 call trsp1(vinvt, xh%lx)
222 type(
cpr_t),
intent(inout) :: cpr
223 real(kind=rp) :: w2(
cpr%Xh%lx,
cpr%Xh%lx,
cpr%Xh%lx)
224 real(kind=rp) :: w1(
cpr%Xh%lx,
cpr%Xh%lx,
cpr%Xh%lx)
225 real(kind=rp) :: specmat(
cpr%Xh%lx,
cpr%Xh%lx)
226 real(kind=rp) :: specmatt(
cpr%Xh%lx,
cpr%Xh%lx)
227 integer :: i, j, k, e, nxyz, nelv
228 character(len=LOG_SIZE) :: log_buf
229 character(len=4) :: space
236 if (
space .eq.
'spec')
then
237 call copy(specmat,
cpr%vinv,
cpr%Xh%lx*
cpr%Xh%lx)
238 call copy(specmatt,
cpr%vinvt,
cpr%Xh%lx*
cpr%Xh%lx)
240 if (
space .eq.
'phys')
then
241 call copy(specmat,
cpr%v,
cpr%Xh%lx*
cpr%Xh%lx)
242 call copy(specmatt,
cpr%vt,
cpr%Xh%lx*
cpr%Xh%lx)
247 if (
space .eq.
'spec')
then
248 call copy(w2,
cpr%fld%x(1,1,1,e), nxyz)
250 call copy(w2,
cpr%fldhat(1,1,1,e), nxyz)
253 call mxm(specmat,
cpr%Xh%lx, w2,
cpr%Xh%lx, w1,
cpr%Xh%lx*
cpr%Xh%lx)
256 call mxm(w1(1,1,k),
cpr%Xh%lx,specmatt,
cpr%Xh%lx,w2(1,1,k),
cpr%Xh%lx)
259 call mxm (w2,
cpr%Xh%lx*
cpr%Xh%lx, specmatt,
cpr%Xh%lx,&
260 cpr%fldhat(1,1,1,e),
cpr%Xh%lx)
267 type(
cpr_t),
intent(inout) :: cpr
268 type(coef_t),
intent(inout) :: coef
269 real(kind=rp) :: w2(
cpr%Xh%lx,
cpr%Xh%lx,
cpr%Xh%lx)
270 real(kind=rp) :: w1(
cpr%Xh%lx,
cpr%Xh%lx,
cpr%Xh%lx)
271 real(kind=rp) :: vsort(
cpr%Xh%lx,
cpr%Xh%lx,
cpr%Xh%lx)
272 real(kind=rp) :: vtrunc(
cpr%Xh%lx,
cpr%Xh%lx,
cpr%Xh%lx)
273 real(kind=rp) :: vtemp(
cpr%Xh%lx,
cpr%Xh%lx,
cpr%Xh%lx)
274 real(kind=rp) :: errvec(
cpr%Xh%lx,
cpr%Xh%lx,
cpr%Xh%lx)
275 real(kind=rp) :: fx(
cpr%Xh%lx,
cpr%Xh%lx)
276 real(kind=rp) :: fy(
cpr%Xh%lx,
cpr%Xh%lx)
277 real(kind=rp) :: fz(
cpr%Xh%lx,
cpr%Xh%lx)
278 real(kind=rp) :: l2norm, oldl2norm, targeterr
279 integer :: isort(cpr%Xh%lx, cpr%Xh%lx, cpr%Xh%lx)
280 integer :: i, j, k, e, nxyz, nelv
281 integer :: kut, kutx, kuty, kutz, nx
282 character(len=LOG_SIZE) :: log_buf
293 call copy(vtemp,
cpr%fldhat(1,1,1,e), nxyz)
294 call copy(vtrunc,
cpr%fldhat(1,1,1,e), nxyz)
301 do while (l2norm .le. targeterr .and. kut .le. (
cpr%Xh%lx-1)*3)
303 call copy(vtrunc, vtemp, nxyz)
313 call copy(w2, vsort, nxyz)
315 call mxm(fx, nx, w2, nx, w1, nx*nx)
318 call mxm(w1(1,1,k), nx, fy, nx, w2(1,1,k), nx)
321 call mxm (w2, nx*nx, fz, nx, vtemp, nx)
325 call reord(vtemp, isort, nxyz)
328 call sub3(errvec,
cpr%fldhat(1,1,1,e), vtemp, nxyz)
336 call copy(
cpr%fldhat(1,1,1,e), vtrunc, nxyz)
343 write(log_buf,
'(A)') &
344 'Debugging info for e=50. Do not forget to delete'
345 call neko_log%message(log_buf)
347 call reord(vsort,isort,nxyz)
349 write(log_buf,
'(A)') &
351 call neko_log%message(log_buf)
353 write(log_buf,
'(A,E15.7)') &
355 call neko_log%message(log_buf)
356 write(log_buf,
'(A)') &
358 call neko_log%message(log_buf)
360 write(log_buf,
'(A,I5)') &
362 call neko_log%message(log_buf)
363 write(log_buf,
'(A)') &
364 'spectral coefficients'
365 call neko_log%message(log_buf)
367 write(log_buf,
'(A,E15.7,A,E15.7)') &
368 'org coeff:', vsort(i,1,1),
' truncated coeff', &
370 call neko_log%message(log_buf)
380 integer,
intent(in) :: nxyz
381 real(kind=rp),
intent(inout) :: vsort(nxyz)
382 real(kind=rp),
intent(inout) :: v(nxyz)
383 integer,
intent(inout) :: isort(nxyz)
384 real(kind=rp) :: wrksort(nxyz)
385 integer :: wrkisort(nxyz)
396 call sort(vsort, isort, nxyz)
399 call flipv(vsort, isort, nxyz)
402 call copy(vsort, v, nxyz)
403 call swap(vsort, isort, nxyz)
409 integer,
intent(in) :: n
410 real(kind=rp),
intent(inout) :: b(n)
411 integer,
intent(inout) :: ind(n)
412 real(kind=rp) :: temp(n)
413 integer :: tempind(n)
430 integer,
intent(in) :: n
431 real(kind=rp),
intent(inout) :: b(n)
432 integer,
intent(inout) :: ind(n)
433 real(kind=rp) :: temp(n)
448 integer,
intent(in) :: n
449 real(kind=rp),
intent(inout) :: b(n)
450 integer,
intent(inout) :: ind(n)
451 real(kind=rp) :: temp(n)
466 integer,
intent(in) :: lx
467 integer,
intent(in) :: kut
468 real(kind=rp),
intent(inout) :: fx(lx,lx)
469 real(kind=rp),
intent(inout) :: fy(lx,lx)
470 real(kind=rp),
intent(inout) :: fz(lx,lx)
471 integer :: i, j, k, k0, kk, kutx, kuty, kutz
479 if (kut .gt. 0 .and. kut .le. (lx-1))
then
484 if (kut .gt. (lx-1) .and. kut .le. (lx-1)*2)
then
489 if (kut .gt. (lx-1)*2 .and. kut .le. (lx-1)*3)
then
531 type(coef_t),
intent(in) :: coef
532 real(kind=rp) :: elemdata(coef%Xh%lx, coef%Xh%lx, coef%Xh%lx)
533 real(kind=rp) :: vole, suma, l2e
534 integer i, e, eg, nxyz
535 character(len=4) ::
space
538 nxyz = coef%Xh%lx*coef%Xh%lx*coef%Xh%lx
541 vole = vole+coef%B(i,1,1,e)
546 if (
space .eq.
'spec')
then
548 suma = suma+elemdata(i,1,1)*elemdata(i,1,1)*coef%jac(i,1,1,e)
551 if (
space .eq.
'phys')
then
553 suma = suma+elemdata(i,1,1)*elemdata(i,1,1)*coef%B(i,1,1,e)
557 l2e = sqrt(suma)/sqrt(vole)
subroutine build_filter_tf(fx, fy, fz, kut, lx)
create filter transfer function
real(kind=rp) function get_elem_l2norm(elemdata, coef, space, e)
subroutine cpr_truncate_wn(cpr, coef)
Truncate the frequencies.
subroutine swap(b, ind, n)
sort the array acording to ind vector
subroutine cpr_init_all(cpr, u)
Initialize cpr.
subroutine sortcoeff(vsort, v, isort, nxyz)
Sort the spectral coefficient in descending order array vsort. The original indices are stored in the...
subroutine flipv(b, ind, n)
Flip vector b and ind.
subroutine cpr_goto_space(cpr, space)
Tranform to spectral space (using tensor product)
subroutine reord(b, ind, n)
reorder the array - inverse of swap
subroutine cpr_generate_specmat(cpr)
Generate spectral tranform matrices.
subroutine, public cpr_free(cpr)
Deallocate coefficients.
Defines a mapping of the degrees of freedom.
subroutine, public copy(a, b, n)
Copy a vector .
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.
Defines a function space.
subroutine, public trsp1(a, n)
In-place transpose of a square tensor.
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
include information needed for compressing fields
The function space for the SEM solution fields.