spectral_error.f90

Go to the documentation of this file.

! Copyright (c) 2022, The Neko Authors
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
!
!   * Redistributions of source code must retain the above copyright
!     notice, this list of conditions and the following disclaimer.
!
!   * Redistributions in binary form must reproduce the above
!     copyright notice, this list of conditions and the following
!     disclaimer in the documentation and/or other materials provided
!     with the distribution.
!
!   * Neither the name of the authors nor the names of its
!     contributors may be used to endorse or promote products derived
!     from this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
module spectral_error
  use num_types, only: rp
  use field, only: field_t
  use coefs, only: coef_t
  use field_list, only: field_list_t
  use math, only: rzero, copy
  use file, only: file_t, file_free
  use time_state, only : time_state_t
  use tensor, only: tnsr3d
  use device_math, only: device_copy
  use gather_scatter
  use neko_config
  use logger, only: neko_log
  use device, only: device_to_host, host_to_device, device_memcpy
  use comm, only: pe_rank
  use utils, only: neko_fname_len, neko_error
  use field_writer, only: field_writer_t
  use simulation_component, only: simulation_component_t
  use json_module, only: json_file
  use json_utils, only: json_get, json_get_or_default
  use case, only: case_t
  use field_registry, only: neko_field_registry
 
  use, intrinsic :: iso_c_binding
  implicit none
  private
 
  type, public, extends(simulation_component_t) :: spectral_error_t
     type(field_t), pointer :: u => null()
     type(field_t), pointer :: v => null()
     type(field_t), pointer :: w => null()
     type(field_t), pointer :: u_hat => null()
     type(field_t), pointer :: v_hat => null()
     type(field_t), pointer :: w_hat => null()
     type(field_t) :: wk
     real(kind=rp) :: seri_small = 1.e-14
     real(kind=rp) :: seri_smallr = 1.e-10
     real(kind=rp) :: seri_smallg = 1.e-5
     real(kind=rp) :: seri_smalls = 0.2
     integer :: seri_np = 4
     integer :: seri_np_max = 4
     integer :: seri_elr = 0
     real(kind=rp), allocatable :: eind_u(:), eind_v(:), eind_w(:)
     real(kind=rp), allocatable :: sig_u(:), sig_v(:), sig_w(:)
     type(field_list_t) :: speri_l
     type(file_t) :: mf_speri
     type(field_writer_t) :: writer
   contains
     procedure, pass(this) :: init => spectral_error_init
     procedure, pass(this) :: free => spectral_error_free
     procedure, pass(this) :: compute_ => spectral_error_compute
     procedure, pass(this) :: get_indicators => spectral_error_get_indicators
 
  end type spectral_error_t
 
contains
 
  subroutine spectral_error_init(this, json, case)
    class(spectral_error_t), intent(inout), target :: this
    type(json_file), intent(inout) :: json
    class(case_t), intent(inout), target :: case
 
    character(len=20) :: fields(3)
 
    ! picks it up. Will also add those fields to the registry.
    fields(1) = "u_hat"
    fields(2) = "v_hat"
    fields(3) = "w_hat"
    call json%add("fields", fields)
 
    call this%init_base(json, case)
    call this%writer%init(json, case)
 
    call spectral_error_init_from_components(this, case%fluid%c_Xh)
 
  end subroutine spectral_error_init
 
  subroutine spectral_error_init_from_components(this, coef)
    class(spectral_error_t), intent(inout) :: this
    type(coef_t), intent(in) :: coef
    integer :: il, jl, aa
    character(len=NEKO_FNAME_LEN) :: fname_speri
 
    this%u => neko_field_registry%get_field("u")
    this%v => neko_field_registry%get_field("v")
    this%w => neko_field_registry%get_field("w")
    this%u_hat => neko_field_registry%get_field("u_hat")
    this%v_hat => neko_field_registry%get_field("v_hat")
    this%w_hat => neko_field_registry%get_field("w_hat")
 
    this%wk = this%u
 
    allocate(this%eind_u(coef%msh%nelv))
    allocate(this%eind_v(coef%msh%nelv))
    allocate(this%eind_w(coef%msh%nelv))
    allocate(this%sig_u(coef%msh%nelv))
    allocate(this%sig_v(coef%msh%nelv))
    allocate(this%sig_w(coef%msh%nelv))
 
    associate(lx1 => coef%Xh%lx, ly1 => coef%Xh%ly, &
         lz1 => coef%Xh%lz, &
         seri_small => this%SERI_SMALL, &
         seri_smallr => this%SERI_SMALLR, &
         seri_smallg => this%SERI_SMALLG, &
         seri_smalls => this%SERI_SMALLS, &
         seri_np => this%SERI_NP, &
         seri_np_max => this%SERI_NP_MAX, &
         seri_elr => this%SERI_ELR &
         )
      ! correctness check
      if (seri_np.gt.seri_np_max) then
         call neko_log%message('SETI_NP greater than SERI_NP_MAX')
      endif
      il = seri_np+seri_elr
      jl = min(lx1,ly1)
      jl = min(jl,lz1)
      if (il.gt.jl) then
         call neko_log%message('SERI_NP+SERI_ELR greater than L?1')
      endif
    end associate
 
  end subroutine spectral_error_init_from_components
 
  subroutine spectral_error_free(this)
    class(spectral_error_t), intent(inout) :: this
 
    if(allocated(this%eind_u)) then
       deallocate(this%eind_u)
    end if
 
    if(allocated(this%eind_v)) then
       deallocate(this%eind_v)
    end if
 
    if(allocated(this%eind_w)) then
       deallocate(this%eind_w)
    end if
 
    if(allocated(this%sig_u)) then
       deallocate(this%sig_u)
    end if
 
    if(allocated(this%sig_w)) then
       deallocate(this%sig_w)
    end if
 
    if(allocated(this%sig_w)) then
       deallocate(this%sig_w)
    end if
 
    call this%wk%free()
 
    nullify(this%u)
    nullify(this%v)
    nullify(this%w)
    nullify(this%u_hat)
    nullify(this%v_hat)
    nullify(this%w_hat)
 
    call this%writer%free()
    call this%free_base()
 
  end subroutine spectral_error_free
 
  subroutine spectral_error_compute(this, time)
    class(spectral_error_t), intent(inout) :: this
    type(time_state_t), intent(in) :: time
 
    integer :: e, i, lx, ly, lz, nelv, n
 
 
    call this%get_indicators(this%case%fluid%c_Xh)
 
    lx = this%u_hat%Xh%lx
    ly = this%u_hat%Xh%ly
    lz = this%u_hat%Xh%lz
    nelv = this%u_hat%msh%nelv
 
    do e = 1,nelv
       do i = 1,lx*ly*lx
          this%u_hat%x(i,1,1,e) = this%eind_u(e)
          this%v_hat%x(i,1,1,e) = this%eind_v(e)
          this%w_hat%x(i,1,1,e) = this%eind_w(e)
       end do
    end do
 
    ! We need this copy to GPU since the field writer does an internal copy
    ! GPU -> CPU before writing the field files. This overwrites the *_hat
    ! host arrays that contain the values.
    if (neko_bcknd_device .eq. 1) then
       call device_memcpy(this%u_hat%x, this%u_hat%x_d, lx*ly*lz*nelv, &
            host_to_device, sync = .true.)
       call device_memcpy(this%v_hat%x, this%v_hat%x_d, lx*ly*lz*nelv, &
            host_to_device, sync = .true.)
       call device_memcpy(this%w_hat%x, this%w_hat%x_d, lx*ly*lz*nelv, &
            host_to_device, sync = .true.)
    end if
 
  end subroutine spectral_error_compute
 
  subroutine transform_to_spec_or_phys(u_hat, u, wk, coef, space)
    type(field_t), intent(inout) :: u_hat
    type(field_t), intent(inout) :: u
    type(field_t), intent(inout) :: wk
    type(coef_t), intent(inout) :: coef
    character(len=4), intent(in) :: space
    integer :: i, j, k, e, nxyz, nelv, n
 
    nxyz = coef%Xh%lx*coef%Xh%lx*coef%Xh%lx
    nelv = coef%msh%nelv
    n = nxyz*nelv
 
    if ((neko_bcknd_hip .eq. 1) .or. (neko_bcknd_cuda .eq. 1) .or. &
         (neko_bcknd_opencl .eq. 1)) then
       call device_copy(wk%x_d, u%x_d, n)
    else
       call copy(wk%x,u%x,n)
    end if
 
    select case(space)
    case('spec')
       call tnsr3d(u_hat%x, coef%Xh%lx, wk%x, &
            coef%Xh%lx,coef%Xh%vinv, &
            coef%Xh%vinvt, coef%Xh%vinvt, nelv)
    case('phys')
       call tnsr3d(u_hat%x, coef%Xh%lx, wk%x, &
            coef%Xh%lx,coef%Xh%v, &
            coef%Xh%vt, coef%Xh%vt, nelv)
    end select
 
    ! Synchronize
    if ((neko_bcknd_hip .eq. 1) .or. (neko_bcknd_cuda .eq. 1) .or. &
         (neko_bcknd_opencl .eq. 1)) then
 
       call device_memcpy(u_hat%x,u_hat%x_d, n, &
            device_to_host, sync=.true.)
    end if
 
  end subroutine transform_to_spec_or_phys
 
  subroutine spectral_error_get_indicators(this,coef)
    class(spectral_error_t), intent(inout) :: this
    type(coef_t), intent(inout) :: coef
    integer :: i
 
    ! Generate the uvwhat field (legendre coeff)
    call transform_to_spec_or_phys(this%u_hat, this%u, this%wk, coef, 'spec')
    call transform_to_spec_or_phys(this%v_hat, this%v, this%wk, coef, 'spec')
    call transform_to_spec_or_phys(this%w_hat, this%w, this%wk, coef, 'spec')
 
    ! Get the spectral error indicator
    call calculate_indicators(this, coef, this%eind_u, this%sig_u, &
         coef%msh%nelv, coef%Xh%lx, coef%Xh%ly, coef%Xh%lz, &
         this%u_hat%x)
    call calculate_indicators(this, coef, this%eind_v, this%sig_v, &
         coef%msh%nelv, coef%Xh%lx, coef%Xh%ly, coef%Xh%lz, &
         this%v_hat%x)
    call calculate_indicators(this, coef, this%eind_w, this%sig_w, &
         coef%msh%nelv, coef%Xh%lx, coef%Xh%ly, coef%Xh%lz, &
         this%w_hat%x)
 
  end subroutine spectral_error_get_indicators
 
  subroutine spectral_error_write(this, t)
    class(spectral_error_t), intent(inout) :: this
    real(kind=rp), intent(in) :: t
 
    integer i, e
    integer lx, ly, lz, nelv
 
    call this%mf_speri%write(this%speri_l,t)
 
  end subroutine spectral_error_write
 
  subroutine calculate_indicators(this, coef, eind, sig, lnelt, LX1, LY1, LZ1, var)
    type(spectral_error_t), intent(inout) :: this
    type(coef_t), intent(in) :: coef
    integer :: lnelt
    integer :: LX1
    integer :: LY1
    integer :: LZ1
    real(kind=rp) :: eind(lnelt)
    real(kind=rp) :: sig(lnelt)
    real(kind=rp) :: var(lx1,ly1,lz1,lnelt)
 
    real(kind=rp) :: xa(coef%Xh%lx,coef%Xh%ly,coef%Xh%lz)
    real(kind=rp) :: xb(coef%Xh%lx,coef%Xh%ly,coef%Xh%lz)
    integer :: i, e
 
    ! zero arrays
    call rzero(eind,lnelt)
    call rzero(sig,lnelt)
 
    ! Get the indicator
    call speri_var(this, eind,sig,var,lnelt,xa,xb, lx1, ly1, lz1)
 
  end subroutine calculate_indicators
 
 
  subroutine speri_var(this, est,sig,var,nell,xa,xb,LX1,LY1,LZ1)
    type(spectral_error_t), intent(inout) :: this
    integer :: nell
    integer :: LX1
    integer :: LY1
    integer :: LZ1
    real(kind=rp) :: est(nell)
    real(kind=rp) :: sig(nell)
    real(kind=rp) :: var(lx1,ly1,lz1,nell)
    real(kind=rp) :: xa(lx1,ly1,lz1)
    real(kind=rp) :: xb(lx1,ly1,lz1)
 
    integer :: il, jl, kl, ll, j_st, j_en, ii
    real(kind=rp) :: coeff(lx1,ly1,lz1)
    real(kind=rp) :: coef11
    real(kind=rp) :: coefx(this%SERI_NP_MAX,ly1,lz1), &
         coefy(this%SERI_NP_MAX,lx1,lz1), &
         coefz(this%SERI_NP_MAX,lx1,ly1)
    real(kind=rp) :: estx, esty, estz
    real(kind=rp) :: sigx, sigy, sigz
    real(kind=rp) :: third
    parameter(third = 1.0/3.0)
 
    do il = 1,nell
       do ii = 1, lx1*ly1*lz1
          coeff(ii,1,1) = var(ii,1,1,il) * var(ii,1,1,il)
       end do
 
       coef11 = coeff(1,1,1)
 
       if (coef11.ge.this%SERI_SMALL) then
          j_st = max(1,lx1-this%SERI_NP+1-this%SERI_ELR)
          j_en = max(1,lx1-this%SERI_ELR)
          do ll = 1,lz1
             do kl = 1,ly1
                do jl = j_st,j_en
                   coefx(j_en-jl+1,kl,ll) = coeff(jl,kl,ll)
                enddo
             enddo
          enddo
          call speri_extrap(this,estx,sigx,coef11,coefx, &
               j_st,j_en,ly1,lz1)
 
          j_st = max(1,ly1-this%SERI_NP+1-this%SERI_ELR)
          j_en = max(1,ly1-this%SERI_ELR)
          do ll = 1,lz1
             do kl = j_st,j_en
                do jl = 1,lx1
                   coefy(j_en-kl+1,jl,ll) = coeff(jl,kl,ll)
                enddo
             enddo
          enddo
          call speri_extrap(this, esty,sigy,coef11,coefy, &
               j_st,j_en,lx1,lz1)
 
          j_st = max(1,lz1-this%SERI_NP+1-this%SERI_ELR)
          j_en = max(1,lz1-this%SERI_ELR)
          do ll = j_st,j_en
             do kl = 1,ly1
                do jl = 1,lx1
                   coefz(j_en-ll+1,jl,kl) = coeff(jl,kl,ll)
                enddo
             enddo
          enddo
          call speri_extrap(this, estz,sigz,coef11,coefz, &
               j_st,j_en,lx1,ly1)
 
          est(il) = sqrt(estx + esty + estz)
          sig(il) = third*(sigx + sigy + sigz)
 
       else
          estx = 0.0
          esty = 0.0
          estz = 0.0
          sigx = -1.0
          sigy = -1.0
          sigz = -1.0
 
          est(il) = 0.0
          sig(il) = -1.0
       endif
 
    end do
 
  end subroutine speri_var
 
  subroutine speri_extrap(this,estx,sigx,coef11,coef, &
       ix_st,ix_en,nyl,nzl)
    implicit none
    type(spectral_error_t), intent(inout) :: this
    integer :: ix_st,ix_en,nyl,nzl
    real(kind=rp) :: coef(this%SERI_NP_MAX,nyl,nzl)
    real(kind=rp) :: coef11
    real(kind=rp) :: estx, sigx
 
    integer :: il, jl, kl, ll ! loop index
    integer :: nsigt, pnr, nzlt
    real(kind=rp) :: sigt, smallr, cmin, cmax, cnm, rtmp, rtmp2, rtmp3
    real(kind=rp) :: sumtmp(4), cffl(this%SERI_NP_MAX)
    real(kind=rp) :: stmp, estt, clog, ctmp, cave, erlog
    logical :: cuse(this%seri_np_max)
 
    associate(seri_small => this%SERI_SMALL, &
         seri_smallr => this%SERI_SMALLR, &
         seri_smallg => this%SERI_SMALLG, &
         seri_smalls => this%SERI_SMALLS, &
         seri_np => this%SERI_NP, &
         seri_np_max => this%SERI_NP_MAX, &
         seri_elr => this%SERI_ELR &
         )
      ! initial values
      estx = 0.0
      sigx = -1.0
 
      ! relative cutoff
      smallr = coef11*seri_smallr
 
      ! number of points
      pnr = ix_en - ix_st +1
 
      ! to few points to interpolate
      !if ((ix_en - ix_st).le.1) return
 
      ! for averaging, initial values
      sigt = 0.0
      nsigt = 0
 
      ! loop over all face points
      nzlt = max(1,nzl - seri_elr) !  for 2D runs
      do il=1,nzlt
         ! weight
         rtmp3 = 1.0/(2.0*(il-1)+1.0)
         do jl=1,nyl - seri_elr
 
            ! find min and max coef along single row
            cffl(1) = coef(1,jl,il)
            cmin = cffl(1)
            cmax = cmin
            do kl =2,pnr
               cffl(kl) = coef(kl,jl,il)
               cmin = min(cmin,cffl(kl))
               cmax = max(cmax,cffl(kl))
            enddo
 
            ! are coefficients sufficiently big
            if((cmin.gt.0.0).and.(cmax.gt.smallr)) then
               ! mark array position we use in iterpolation
               do kl =1,pnr
                  cuse(kl) = .true.
               enddo
               ! max n for polynomial order
               cnm = real(ix_en)
 
               ! check if all the points should be taken into account
               ! in original code by Catherine Mavriplis this part is written
               ! for 4 points, so I place if statement first
               if (pnr.eq.4) then
                  ! should we neglect last values
                  if ((cffl(1).lt.smallr).and. &
                       (cffl(2).lt.smallr)) then
                     if (cffl(3).lt.smallr) then
                        cuse(1) = .false.
                        cuse(2) = .false.
                        cnm = real(ix_en-2)
                     else
                        cuse(1) = .false.
                        cnm = real(ix_en-1)
                     endif
                  else
                     ! should we take stronger gradient
                     if ((cffl(1)/cffl(2).lt.seri_smallg).and. &
                          (cffl(3)/cffl(4).lt.seri_smallg)) then
                        cuse(1) = .false.
                        cuse(3) = .false.
                        cnm = real(ix_en-1)
                     elseif ((cffl(2)/cffl(1).lt.seri_smallg).and. &
                          (cffl(4)/cffl(3).lt.seri_smallg)) then
                        cuse(2) = .false.
                        cuse(4) = .false.
                     endif
                  endif
               endif
 
               ! get sigma for given face point
               do kl =1,4
                  sumtmp(kl) = 0.0
               enddo
               ! find new min and count number of points
               cmin = cmax
               cmax = 0.0
               do kl =1,pnr
                  if(cuse(kl)) then
                     rtmp = real(ix_en-kl)
                     rtmp2 = log(cffl(kl))
                     sumtmp(1) = sumtmp(1) +rtmp2
                     sumtmp(2) = sumtmp(2) +rtmp
                     sumtmp(3) = sumtmp(3) +rtmp*rtmp
                     sumtmp(4) = sumtmp(4) +rtmp2*rtmp
                     ! find new min and count used points
                     cmin = min(cffl(kl),cmin)
                     cmax = cmax + 1.0
                  endif
               enddo
               ! decay rate along single row
               stmp = (sumtmp(1)*sumtmp(2) - sumtmp(4)*cmax)/ &
                    (sumtmp(3)*cmax - sumtmp(2)*sumtmp(2))
               ! for averaging
               sigt = sigt + stmp
               nsigt = nsigt + 1
 
               ! get error estimator depending on calculated decay rate
               estt = 0.0
               if (stmp.lt.seri_smalls) then
                  estt = cmin
               else
                  ! get averaged constant in front of c*exp(-sig*n)
                  clog = (sumtmp(1)+stmp*sumtmp(2))/cmax
                  ctmp = exp(clog)
                  ! average exponent
                  cave = sumtmp(1)/cmax
                  ! check quality of approximation comparing is to the constant cave
                  do kl =1,2
                     sumtmp(kl) = 0.0
                  enddo
                  do kl =1,pnr
                     if(cuse(kl)) then
                        erlog = clog - stmp*real(ix_en-kl)
                        sumtmp(1) = sumtmp(1)+ &
                             (erlog-log(cffl(kl)))**2
                        sumtmp(2) = sumtmp(2)+ &
                             (erlog-cave)**2
                     endif
                  enddo
                  rtmp = 1.0 - sumtmp(1)/sumtmp(2)
                  if (rtmp.lt.seri_smalls) then
                     estt = cmin
                  else
                     ! last coefficient is not included in error estimator
                     estt = ctmp/stmp*exp(-stmp*cnm)
                  endif
               endif
               ! add contribution to error estimator; variable weight
               estx = estx + estt/(2.0*(jl-1)+1.0)*rtmp3
            endif ! if((cmin.gt.0.0).and.(cmax.gt.smallr))
         enddo
      enddo
      ! constant weight
      ! Multiplication by 4 in 2D / 8 in 3D
      ! Normalization of the error by the volume of the reference element
      ! which is equal to 4 in 2D / 8 in 3D
      ! ==> Both operations cancel each other
      estx = estx/(2.0*(ix_en-1)+1.0)
 
      ! final everaging
      ! sigt = 2*sigma so we divide by 2
      if (nsigt.gt.0) then
         sigx = 0.5*sigt/nsigt
      endif
 
    end associate
 
  end subroutine speri_extrap
 
end module spectral_error