48 use,
intrinsic :: iso_c_binding
69 real(kind=
rp) :: seri_small = 1.e-14
71 real(kind=
rp) :: seri_smallr = 1.e-10
73 real(kind=
rp) :: seri_smallg = 1.e-5
75 real(kind=
rp) :: seri_smalls = 0.2
77 integer :: seri_np = 4
78 integer :: seri_np_max = 4
80 integer :: seri_elr = 0
82 real(kind=
rp),
allocatable :: eind_u(:), eind_v(:), eind_w(:)
84 real(kind=
rp),
allocatable :: sig_u(:), sig_v(:), sig_w(:)
110 type(
field_t),
intent(in),
target :: u
111 type(
field_t),
intent(in),
target :: v
112 type(
field_t),
intent(in),
target :: w
113 type(
coef_t),
intent(in) :: coef
114 integer :: il, jl, aa
115 character(len=NEKO_FNAME_LEN) :: fname_speri
130 allocate(this%eind_u(coef%msh%nelv))
131 allocate(this%eind_v(coef%msh%nelv))
132 allocate(this%eind_w(coef%msh%nelv))
133 allocate(this%sig_u(coef%msh%nelv))
134 allocate(this%sig_v(coef%msh%nelv))
135 allocate(this%sig_w(coef%msh%nelv))
138 associate(lx1 => coef%Xh%lx, ly1 => coef%Xh%ly, &
140 seri_small => this%SERI_SMALL, &
141 seri_smallr => this%SERI_SMALLR, &
142 seri_smallg => this%SERI_SMALLG, &
143 seri_smalls => this%SERI_SMALLS, &
144 seri_np => this%SERI_NP, &
145 seri_np_max => this%SERI_NP_MAX, &
146 seri_elr => this%SERI_ELR &
149 if (seri_np.gt.seri_np_max)
then
150 if (
pe_rank.eq.0)
write(*,*)
'SETI_NP greater than SERI_NP_MAX'
152 il = seri_np+seri_elr
156 if (
pe_rank.eq.0)
write(*,*)
'SERI_NP+SERI_ELR greater than L?1'
161 call list_init3(this%speri_l,this%u_hat,this%v_hat, &
164 fname_speri =
'speri.fld'
165 this%mf_speri =
file_t(fname_speri)
174 if(
allocated(this%eind_u))
then
175 deallocate(this%eind_u)
178 if(
allocated(this%eind_v))
then
179 deallocate(this%eind_v)
182 if(
allocated(this%eind_w))
then
183 deallocate(this%eind_w)
186 if(
allocated(this%sig_u))
then
187 deallocate(this%sig_u)
190 if(
allocated(this%sig_w))
then
191 deallocate(this%sig_w)
194 if(
allocated(this%sig_w))
then
195 deallocate(this%sig_w)
198 call this%u_hat%free()
199 call this%v_hat%free()
200 call this%w_hat%free()
209 call file_free(this%mf_speri)
221 type(field_t),
intent(inout) :: u_hat
222 type(field_t),
intent(inout) :: u
223 type(field_t),
intent(inout) :: wk
224 type(coef_t),
intent(inout) :: coef
225 character(len=4),
intent(in) :: space
226 integer :: i, j, k, e, nxyz, nelv, n
229 nxyz = coef%Xh%lx*coef%Xh%lx*coef%Xh%lx
234 if ((neko_bcknd_hip .eq. 1) .or. (neko_bcknd_cuda .eq. 1) .or. &
235 (neko_bcknd_opencl .eq. 1))
then
236 call device_copy(wk%x_d, u%x_d, n)
238 call copy(wk%x,u%x,n)
243 call tnsr3d(u_hat%x, coef%Xh%lx, wk%x, &
244 coef%Xh%lx,coef%Xh%vinv, &
245 coef%Xh%vinvt, coef%Xh%vinvt, nelv)
247 call tnsr3d(u_hat%x, coef%Xh%lx, wk%x, &
248 coef%Xh%lx,coef%Xh%v, &
249 coef%Xh%vt, coef%Xh%vt, nelv)
253 if ((neko_bcknd_hip .eq. 1) .or. (neko_bcknd_cuda .eq. 1) .or. &
254 (neko_bcknd_opencl .eq. 1))
then
256 call device_memcpy(u_hat%x,u_hat%x_d, n, &
257 device_to_host, sync=.true.)
266 type(coef_t),
intent(inout) :: coef
277 coef%Xh%lx, coef%Xh%ly, coef%Xh%lz, &
280 coef%Xh%lx, coef%Xh%ly, coef%Xh%lz, &
283 coef%Xh%lx, coef%Xh%ly, coef%Xh%lz, &
292 real(kind=rp),
intent(in) :: t
295 integer lx, ly, lz, nelv
297 lx = this%u_hat%Xh%lx
298 ly = this%u_hat%Xh%ly
299 lz = this%u_hat%Xh%lz
300 nelv = this%u_hat%msh%nelv
305 this%u_hat%x(i,1,1,e) = this%eind_u(e)
306 this%v_hat%x(i,1,1,e) = this%eind_v(e)
307 this%w_hat%x(i,1,1,e) = this%eind_w(e)
314 call this%mf_speri%write(this%speri_l,t)
329 type(coef_t),
intent(in) :: coef
334 real(kind=rp) :: eind(lnelt)
335 real(kind=rp) :: sig(lnelt)
336 real(kind=rp) :: var(lx1,ly1,lz1,lnelt)
338 real(kind=rp) :: xa(coef%Xh%lx,coef%Xh%ly,coef%Xh%lz)
339 real(kind=rp) :: xb(coef%Xh%lx,coef%Xh%ly,coef%Xh%lz)
343 call rzero(eind,lnelt)
344 call rzero(sig,lnelt)
347 call speri_var(this, eind,sig,var,lnelt,xa,xb, lx1, ly1, lz1)
362 subroutine speri_var(this, est,sig,var,nell,xa,xb,LX1,LY1,LZ1)
368 real(kind=rp) :: est(nell)
369 real(kind=rp) :: sig(nell)
370 real(kind=rp) :: var(lx1,ly1,lz1,nell)
371 real(kind=rp) :: xa(lx1,ly1,lz1)
372 real(kind=rp) :: xb(lx1,ly1,lz1)
375 integer :: il, jl, kl, ll, j_st, j_en, ii
377 real(kind=rp) :: coeff(lx1,ly1,lz1)
379 real(kind=rp) :: coef11
381 real(kind=rp) :: coefx(this%SERI_NP_MAX,ly1,lz1), &
382 coefy(this%SERI_NP_MAX,lx1,lz1), &
383 coefz(this%SERI_NP_MAX,lx1,ly1)
385 real(kind=rp) :: estx, esty, estz
387 real(kind=rp) :: sigx, sigy, sigz
388 real(kind=rp) :: third
389 parameter(third = 1.0/3.0)
395 do ii = 1, lx1*ly1*lz1
396 coeff(ii,1,1) = var(ii,1,1,il) * var(ii,1,1,il)
400 coef11 = coeff(1,1,1)
403 if (coef11.ge.this%SERI_SMALL)
then
408 j_st =
max(1,lx1-this%SERI_NP+1-this%SERI_ELR)
409 j_en =
max(1,lx1-this%SERI_ELR)
413 coefx(j_en-jl+1,kl,ll) = coeff(jl,kl,ll)
424 j_st =
max(1,ly1-this%SERI_NP+1-this%SERI_ELR)
425 j_en =
max(1,ly1-this%SERI_ELR)
429 coefy(j_en-kl+1,jl,ll) = coeff(jl,kl,ll)
440 j_st =
max(1,lz1-this%SERI_NP+1-this%SERI_ELR)
441 j_en =
max(1,lz1-this%SERI_ELR)
445 coefz(j_en-ll+1,jl,kl) = coeff(jl,kl,ll)
454 est(il) = sqrt(estx + esty + estz)
455 sig(il) = third*(sigx + sigy + sigz)
490 integer :: ix_st,ix_en,nyl,nzl
492 real(kind=rp) :: coef(this%SERI_NP_MAX,nyl,nzl)
494 real(kind=rp) :: coef11
496 real(kind=rp) :: estx, sigx
499 integer :: il, jl, kl, ll
500 integer :: nsigt, pnr, nzlt
501 real(kind=rp) :: sigt, smallr, cmin, cmax, cnm, rtmp, rtmp2, rtmp3
502 real(kind=rp) :: sumtmp(4), cffl(this%SERI_NP_MAX)
503 real(kind=rp) :: stmp, estt, clog, ctmp, cave, erlog
504 logical :: cuse(this%seri_np_max)
506 associate(seri_small => this%SERI_SMALL, &
507 seri_smallr => this%SERI_SMALLR, &
508 seri_smallg => this%SERI_SMALLG, &
509 seri_smalls => this%SERI_SMALLS, &
510 seri_np => this%SERI_NP, &
511 seri_np_max => this%SERI_NP_MAX, &
512 seri_elr => this%SERI_ELR &
519 smallr = coef11*seri_smallr
522 pnr = ix_en - ix_st +1
532 nzlt =
max(1,nzl - seri_elr)
535 rtmp3 = 1.0/(2.0*(il-1)+1.0)
536 do jl=1,nyl - seri_elr
539 cffl(1) = coef(1,jl,il)
543 cffl(kl) = coef(kl,jl,il)
544 cmin = min(cmin,cffl(kl))
545 cmax =
max(cmax,cffl(kl))
549 if((cmin.gt.0.0).and.(cmax.gt.smallr))
then
562 if ((cffl(1).lt.smallr).and. &
563 (cffl(2).lt.smallr))
then
564 if (cffl(3).lt.smallr)
then
574 if ((cffl(1)/cffl(2).lt.seri_smallg).and. &
575 (cffl(3)/cffl(4).lt.seri_smallg))
then
579 elseif ((cffl(2)/cffl(1).lt.seri_smallg).and. &
580 (cffl(4)/cffl(3).lt.seri_smallg))
then
596 rtmp =
real(ix_en-kl)
597 rtmp2 = log(cffl(kl))
598 sumtmp(1) = sumtmp(1) +rtmp2
599 sumtmp(2) = sumtmp(2) +rtmp
600 sumtmp(3) = sumtmp(3) +rtmp*rtmp
601 sumtmp(4) = sumtmp(4) +rtmp2*rtmp
603 cmin = min(cffl(kl),cmin)
608 stmp = (sumtmp(1)*sumtmp(2) - sumtmp(4)*cmax)/ &
609 (sumtmp(3)*cmax - sumtmp(2)*sumtmp(2))
616 if (stmp.lt.seri_smalls)
then
620 clog = (sumtmp(1)+stmp*sumtmp(2))/cmax
623 cave = sumtmp(1)/cmax
630 erlog = clog - stmp*
real(ix_en-kl)
631 sumtmp(1) = sumtmp(1)+ &
632 (erlog-log(cffl(kl)))**2
633 sumtmp(2) = sumtmp(2)+ &
637 rtmp = 1.0 - sumtmp(1)/sumtmp(2)
638 if (rtmp.lt.seri_smalls)
then
642 estt = ctmp/stmp*exp(-stmp*cnm)
646 estx = estx + estt/(2.0*(jl-1)+1.0)*rtmp3
655 estx = estx/(2.0*(ix_en-1)+1.0)
660 sigx = 0.5*sigt/nsigt
673 type(field_list_t),
intent(inout) :: list
674 type(field_t) ,
target:: uu
675 type(field_t) ,
target:: vv
676 type(field_t) ,
target:: ww
678 allocate(list%fields(3))
679 list%fields(1)%f => uu
680 list%fields(2)%f => vv
681 list%fields(3)%f => ww
688 type(field_list_t),
intent(inout) :: list
690 deallocate(list%fields)
Copy data between host and device (or device and device)
subroutine, public device_copy(a_d, b_d, n)
Device abstraction, common interface for various accelerators.
integer, parameter, public device_to_host
Module for file I/O operations.
subroutine file_free(this)
File operation destructor.
subroutine, public copy(a, b, n)
Copy a vector .
subroutine, public rzero(a, n)
Zero a real vector.
integer, parameter, public rp
Global precision used in computations.
Defines a function space.
Implements type spectral_error_indicator_t.
subroutine speri_extrap(this, estx, sigx, coef11, coef, ix_st, ix_en, nyl, nzl)
Extrapolate the Legendre spectrum from the last points.
subroutine transform_to_spec_or_phys(u_hat, u, wk, coef, space)
Transform a field u > u_hat into physical or spectral space the result of the transformation is in u_...
subroutine list_final3(list)
Helper function to finalize field list to write.
subroutine spec_err_ind_init(this, u, v, w, coef)
Constructor.
subroutine spec_err_ind_get(this, coef)
Transform and get the spectral error indicators.
subroutine list_init3(list, uu, vv, ww)
Helper function to initialize field list to write.
subroutine spec_err_ind_write(this, t)
Write error indicators in a field file.
subroutine speri_var(this, est, sig, var, nell, xa, xb, LX1, LY1, LZ1)
Calculate the indicator in a specified variable.
subroutine calculate_indicators(this, coef, eind, sig, lnelt, LX1, LY1, LZ1, var)
Wrapper for old fortran 77 subroutines.
subroutine spec_err_ind_free(this)
Destructor.
subroutine, public tnsr3d(v, nv, u, nu, A, Bt, Ct, nelv)
Tensor product performed on nelv elements.
integer, parameter neko_fname_len
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
field_list_t, To be able to group fields together
A wrapper around a polymorphic generic_file_t that handles its init. This is essentially a factory fo...
Provides tools to calculate the spectral error indicator.