51 use json_module,
only: json_file
56 use,
intrinsic :: iso_c_binding
77 real(kind=
rp) :: seri_small = 1.e-14
79 real(kind=
rp) :: seri_smallr = 1.e-10
81 real(kind=
rp) :: seri_smallg = 1.e-5
83 real(kind=
rp) :: seri_smalls = 0.2
85 integer :: seri_np = 4
86 integer :: seri_np_max = 4
88 integer :: seri_elr = 0
90 real(kind=
rp),
allocatable :: eind_u(:), eind_v(:), eind_w(:)
92 real(kind=
rp),
allocatable :: sig_u(:), sig_v(:), sig_w(:)
116 type(json_file),
intent(inout) :: json
117 class(
case_t),
intent(inout),
target :: case
119 character(len=20) :: fields(3)
126 call json%add(
"fields", fields)
128 call this%init_base(json,
case)
129 call this%writer%init(json,
case)
139 type(
coef_t),
intent(in) :: coef
140 integer :: il, jl, aa
141 character(len=NEKO_FNAME_LEN) :: fname_speri
154 allocate(this%eind_u(coef%msh%nelv))
155 allocate(this%eind_v(coef%msh%nelv))
156 allocate(this%eind_w(coef%msh%nelv))
157 allocate(this%sig_u(coef%msh%nelv))
158 allocate(this%sig_v(coef%msh%nelv))
159 allocate(this%sig_w(coef%msh%nelv))
162 associate(lx1 => coef%Xh%lx, ly1 => coef%Xh%ly, &
164 seri_small => this%SERI_SMALL, &
165 seri_smallr => this%SERI_SMALLR, &
166 seri_smallg => this%SERI_SMALLG, &
167 seri_smalls => this%SERI_SMALLS, &
168 seri_np => this%SERI_NP, &
169 seri_np_max => this%SERI_NP_MAX, &
170 seri_elr => this%SERI_ELR &
173 if (seri_np.gt.seri_np_max)
then
174 call neko_log%message(
'SETI_NP greater than SERI_NP_MAX')
176 il = seri_np+seri_elr
180 call neko_log%message(
'SERI_NP+SERI_ELR greater than L?1')
190 if(
allocated(this%eind_u))
then
191 deallocate(this%eind_u)
194 if(
allocated(this%eind_v))
then
195 deallocate(this%eind_v)
198 if(
allocated(this%eind_w))
then
199 deallocate(this%eind_w)
202 if(
allocated(this%sig_u))
then
203 deallocate(this%sig_u)
206 if(
allocated(this%sig_w))
then
207 deallocate(this%sig_w)
210 if(
allocated(this%sig_w))
then
211 deallocate(this%sig_w)
223 call this%writer%free()
224 call this%free_base()
231 real(kind=rp),
intent(in) :: t
232 integer,
intent(in) :: tstep
234 integer :: e, i, lx, ly, lz, nelv, n
237 call this%get_indicators(this%case%fluid%c_Xh)
239 lx = this%u_hat%Xh%lx
240 ly = this%u_hat%Xh%ly
241 lz = this%u_hat%Xh%lz
242 nelv = this%u_hat%msh%nelv
247 this%u_hat%x(i,1,1,e) = this%eind_u(e)
248 this%v_hat%x(i,1,1,e) = this%eind_v(e)
249 this%w_hat%x(i,1,1,e) = this%eind_w(e)
256 if (neko_bcknd_device .eq. 1)
then
257 call device_memcpy(this%u_hat%x, this%u_hat%x_d, lx*ly*lz*nelv, &
258 host_to_device, sync = .true.)
259 call device_memcpy(this%v_hat%x, this%v_hat%x_d, lx*ly*lz*nelv, &
260 host_to_device, sync = .true.)
261 call device_memcpy(this%w_hat%x, this%w_hat%x_d, lx*ly*lz*nelv, &
262 host_to_device, sync = .true.)
275 type(field_t),
intent(inout) :: u_hat
276 type(field_t),
intent(inout) :: u
277 type(field_t),
intent(inout) :: wk
278 type(coef_t),
intent(inout) :: coef
279 character(len=4),
intent(in) :: space
280 integer :: i, j, k, e, nxyz, nelv, n
283 nxyz = coef%Xh%lx*coef%Xh%lx*coef%Xh%lx
288 if ((neko_bcknd_hip .eq. 1) .or. (neko_bcknd_cuda .eq. 1) .or. &
289 (neko_bcknd_opencl .eq. 1))
then
290 call device_copy(wk%x_d, u%x_d, n)
292 call copy(wk%x,u%x,n)
297 call tnsr3d(u_hat%x, coef%Xh%lx, wk%x, &
298 coef%Xh%lx,coef%Xh%vinv, &
299 coef%Xh%vinvt, coef%Xh%vinvt, nelv)
301 call tnsr3d(u_hat%x, coef%Xh%lx, wk%x, &
302 coef%Xh%lx,coef%Xh%v, &
303 coef%Xh%vt, coef%Xh%vt, nelv)
307 if ((neko_bcknd_hip .eq. 1) .or. (neko_bcknd_cuda .eq. 1) .or. &
308 (neko_bcknd_opencl .eq. 1))
then
310 call device_memcpy(u_hat%x,u_hat%x_d, n, &
311 device_to_host, sync=.true.)
320 type(coef_t),
intent(inout) :: coef
330 coef%msh%nelv, coef%Xh%lx, coef%Xh%ly, coef%Xh%lz, &
333 coef%msh%nelv, coef%Xh%lx, coef%Xh%ly, coef%Xh%lz, &
336 coef%msh%nelv, coef%Xh%lx, coef%Xh%ly, coef%Xh%lz, &
345 real(kind=rp),
intent(in) :: t
348 integer lx, ly, lz, nelv
353 call this%mf_speri%write(this%speri_l,t)
368 type(coef_t),
intent(in) :: coef
373 real(kind=rp) :: eind(lnelt)
374 real(kind=rp) :: sig(lnelt)
375 real(kind=rp) :: var(lx1,ly1,lz1,lnelt)
377 real(kind=rp) :: xa(coef%Xh%lx,coef%Xh%ly,coef%Xh%lz)
378 real(kind=rp) :: xb(coef%Xh%lx,coef%Xh%ly,coef%Xh%lz)
382 call rzero(eind,lnelt)
383 call rzero(sig,lnelt)
386 call speri_var(this, eind,sig,var,lnelt,xa,xb, lx1, ly1, lz1)
401 subroutine speri_var(this, est,sig,var,nell,xa,xb,LX1,LY1,LZ1)
407 real(kind=rp) :: est(nell)
408 real(kind=rp) :: sig(nell)
409 real(kind=rp) :: var(lx1,ly1,lz1,nell)
410 real(kind=rp) :: xa(lx1,ly1,lz1)
411 real(kind=rp) :: xb(lx1,ly1,lz1)
414 integer :: il, jl, kl, ll, j_st, j_en, ii
416 real(kind=rp) :: coeff(lx1,ly1,lz1)
418 real(kind=rp) :: coef11
420 real(kind=rp) :: coefx(this%SERI_NP_MAX,ly1,lz1), &
421 coefy(this%SERI_NP_MAX,lx1,lz1), &
422 coefz(this%SERI_NP_MAX,lx1,ly1)
424 real(kind=rp) :: estx, esty, estz
426 real(kind=rp) :: sigx, sigy, sigz
427 real(kind=rp) :: third
428 parameter(third = 1.0/3.0)
434 do ii = 1, lx1*ly1*lz1
435 coeff(ii,1,1) = var(ii,1,1,il) * var(ii,1,1,il)
439 coef11 = coeff(1,1,1)
442 if (coef11.ge.this%SERI_SMALL)
then
447 j_st =
max(1,lx1-this%SERI_NP+1-this%SERI_ELR)
448 j_en =
max(1,lx1-this%SERI_ELR)
452 coefx(j_en-jl+1,kl,ll) = coeff(jl,kl,ll)
463 j_st =
max(1,ly1-this%SERI_NP+1-this%SERI_ELR)
464 j_en =
max(1,ly1-this%SERI_ELR)
468 coefy(j_en-kl+1,jl,ll) = coeff(jl,kl,ll)
479 j_st =
max(1,lz1-this%SERI_NP+1-this%SERI_ELR)
480 j_en =
max(1,lz1-this%SERI_ELR)
484 coefz(j_en-ll+1,jl,kl) = coeff(jl,kl,ll)
493 est(il) = sqrt(estx + esty + estz)
494 sig(il) = third*(sigx + sigy + sigz)
528 integer :: ix_st,ix_en,nyl,nzl
530 real(kind=rp) :: coef(this%SERI_NP_MAX,nyl,nzl)
532 real(kind=rp) :: coef11
534 real(kind=rp) :: estx, sigx
537 integer :: il, jl, kl, ll
538 integer :: nsigt, pnr, nzlt
539 real(kind=rp) :: sigt, smallr, cmin, cmax, cnm, rtmp, rtmp2, rtmp3
540 real(kind=rp) :: sumtmp(4), cffl(this%SERI_NP_MAX)
541 real(kind=rp) :: stmp, estt, clog, ctmp, cave, erlog
542 logical :: cuse(this%seri_np_max)
544 associate(seri_small => this%SERI_SMALL, &
545 seri_smallr => this%SERI_SMALLR, &
546 seri_smallg => this%SERI_SMALLG, &
547 seri_smalls => this%SERI_SMALLS, &
548 seri_np => this%SERI_NP, &
549 seri_np_max => this%SERI_NP_MAX, &
550 seri_elr => this%SERI_ELR &
557 smallr = coef11*seri_smallr
560 pnr = ix_en - ix_st +1
570 nzlt =
max(1,nzl - seri_elr)
573 rtmp3 = 1.0/(2.0*(il-1)+1.0)
574 do jl=1,nyl - seri_elr
577 cffl(1) = coef(1,jl,il)
581 cffl(kl) = coef(kl,jl,il)
582 cmin = min(cmin,cffl(kl))
583 cmax =
max(cmax,cffl(kl))
587 if((cmin.gt.0.0).and.(cmax.gt.smallr))
then
600 if ((cffl(1).lt.smallr).and. &
601 (cffl(2).lt.smallr))
then
602 if (cffl(3).lt.smallr)
then
612 if ((cffl(1)/cffl(2).lt.seri_smallg).and. &
613 (cffl(3)/cffl(4).lt.seri_smallg))
then
617 elseif ((cffl(2)/cffl(1).lt.seri_smallg).and. &
618 (cffl(4)/cffl(3).lt.seri_smallg))
then
634 rtmp =
real(ix_en-kl)
635 rtmp2 = log(cffl(kl))
636 sumtmp(1) = sumtmp(1) +rtmp2
637 sumtmp(2) = sumtmp(2) +rtmp
638 sumtmp(3) = sumtmp(3) +rtmp*rtmp
639 sumtmp(4) = sumtmp(4) +rtmp2*rtmp
641 cmin = min(cffl(kl),cmin)
646 stmp = (sumtmp(1)*sumtmp(2) - sumtmp(4)*cmax)/ &
647 (sumtmp(3)*cmax - sumtmp(2)*sumtmp(2))
654 if (stmp.lt.seri_smalls)
then
658 clog = (sumtmp(1)+stmp*sumtmp(2))/cmax
661 cave = sumtmp(1)/cmax
668 erlog = clog - stmp*
real(ix_en-kl)
669 sumtmp(1) = sumtmp(1)+ &
670 (erlog-log(cffl(kl)))**2
671 sumtmp(2) = sumtmp(2)+ &
675 rtmp = 1.0 - sumtmp(1)/sumtmp(2)
676 if (rtmp.lt.seri_smalls)
then
680 estt = ctmp/stmp*exp(-stmp*cnm)
684 estx = estx + estt/(2.0*(jl-1)+1.0)*rtmp3
693 estx = estx/(2.0*(ix_en-1)+1.0)
698 sigx = 0.5*sigt/nsigt
Copy data between host and device (or device and device)
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.
Defines a simulation case.
subroutine, public device_copy(a_d, b_d, n)
Copy a vector .
Device abstraction, common interface for various accelerators.
integer, parameter, public host_to_device
integer, parameter, public device_to_host
Defines a registry for storing solution fields.
type(field_registry_t), target, public neko_field_registry
Global field registry.
Implements the field_writer_t type.
Module for file I/O operations.
subroutine file_free(this)
File operation destructor.
Utilities for retrieving parameters from the case files.
type(log_t), public neko_log
Global log stream.
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.
Simulation components are objects that encapsulate functionality that can be fit to a particular comp...
Defines a function space.
Implements type spectral_error_t.
subroutine speri_extrap(this, estx, sigx, coef11, coef, ix_st, ix_en, nyl, nzl)
Extrapolate the Legendre spectrum from the last points.
subroutine spectral_error_write(this, t)
Write error indicators in a field file.
subroutine spectral_error_free(this)
Destructor.
subroutine spectral_error_get_indicators(this, coef)
Transform and get the spectral error indicators.
subroutine spectral_error_init(this, json, case)
Constructor.
subroutine spectral_error_compute(this, t, tstep)
Compute the spectral error indicator.
subroutine transform_to_spec_or_phys(u_hat, u, wk, coef, space)
Transform a field u to u_hat into physical or spectral space the result of the transformation is in u...
subroutine calculate_indicators(this, coef, eind, sig, lnelt, LX1, LY1, LZ1, var)
Wrapper for old fortran 77 subroutines.
subroutine spectral_error_init_from_attributes(this, coef)
Actual constructor.
subroutine speri_var(this, est, sig, var, nell, xa, xb, LX1, LY1, LZ1)
Calculate the indicator in a specified variable.
subroutine, public tnsr3d(v, nv, u, nu, A, Bt, Ct, nelv)
Tensor product performed on nelv elements.
integer, parameter, public 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 simulation component that writes a 3d field to a file.
A wrapper around a polymorphic generic_file_t that handles its init. This is essentially a factory fo...
Base abstract class for simulation components.
Provides tools to calculate the spectral error indicator.