Estimate_RefinementIndicator2D_t Subroutine

  subroutine Estimate_RefinementIndicator2D_t(this,solution,ivar)
    !! Compute the per-element modal-energy indicator sigma_e and refine/keep/coarsen flag from
    !! the nodal solution field. ivar selects the driving variable in [1,solution%nVar]; passing
    !! SELF_AMR_ALLVARS (=0) reduces the indicator over all variables by taking, per element, the
    !! largest (least smooth) smoothness ratio S_e before the log10.
    implicit none
    class(RefinementIndicator2D_t),intent(inout) :: this
    class(Scalar2D),intent(in) :: solution
    integer,intent(in) :: ivar
    ! Local
    integer :: iel
    real(prec),parameter :: energyFloor = epsilon(1.0_prec)

    if(solution%N /= this%N) then
      print*,__FILE__,':',__LINE__, &
        ' : Error : solution degree does not match indicator degree.'
      stop 1
    endif
    if(solution%nElem /= this%nElem) then
      print*,__FILE__,':',__LINE__, &
        ' : Error : solution element count does not match indicator element count.'
      stop 1
    endif
    if(ivar < 0 .or. ivar > solution%nVar) then
      print*,__FILE__,':',__LINE__, &
        ' : Error : driving-variable index out of range.'
      stop 1
    endif

    do concurrent(iel=1:this%nElem)
      block
        integer :: j,p,q,ii,v,v0,v1,Np
        real(prec) :: tmp(1:this%N+1,1:this%N+1)
        real(prec) :: uhat(1:this%N+1,1:this%N+1)
        real(prec) :: acc,etot,eclip1,eclip2,r1,r2,se,semax

        Np = this%N+1
        if(ivar == SELF_AMR_ALLVARS) then
          v0 = 1
          v1 = solution%nVar
        else
          v0 = ivar
          v1 = ivar
        endif

        semax = 0.0_prec
        do v = v0,v1
          ! Pass 1 (xi direction): tmp(p,j) = sum_i Pmodal(i,p) * u(i,j)
          do j = 1,Np
            do p = 1,Np
              acc = 0.0_prec
              do ii = 1,Np
                acc = acc+this%Pmodal(ii,p)*solution%interior(ii,j,iel,v)
              enddo
              tmp(p,j) = acc
            enddo
          enddo
          ! Pass 2 (eta direction): uhat(p,q) = sum_j Pmodal(j,q) * tmp(p,j)
          do q = 1,Np
            do p = 1,Np
              acc = 0.0_prec
              do ii = 1,Np
                acc = acc+this%Pmodal(ii,q)*tmp(p,ii)
              enddo
              uhat(p,q) = acc
            enddo
          enddo

          ! Modal energies. eclip1 drops the highest mode in each direction,
          ! eclip2 drops the highest two (only meaningful for N >= 2).
          etot = 0.0_prec
          eclip1 = 0.0_prec
          eclip2 = 0.0_prec
          do q = 1,Np
            do p = 1,Np
              etot = etot+uhat(p,q)*uhat(p,q)
              if(p <= Np-1 .and. q <= Np-1) eclip1 = eclip1+uhat(p,q)*uhat(p,q)
              if(p <= Np-2 .and. q <= Np-2) eclip2 = eclip2+uhat(p,q)*uhat(p,q)
            enddo
          enddo

          if(etot <= energyFloor) then
            ! Field is (near) identically zero on this element: perfectly resolved.
            se = 0.0_prec
          else
            r1 = (etot-eclip1)/etot
            if(this%N >= 2 .and. eclip1 > energyFloor) then
              r2 = (eclip1-eclip2)/eclip1
            else
              r2 = 0.0_prec
            endif
            se = max(r1,r2)
          endif

          semax = max(semax,se)
        enddo

        this%indicator(iel) = log10(max(semax,energyFloor))
        if(this%indicator(iel) > this%refineThreshold) then
          this%flag(iel) = SELF_AMR_REFINE
        elseif(this%indicator(iel) < this%coarsenThreshold) then
          this%flag(iel) = SELF_AMR_COARSEN
        else
          this%flag(iel) = SELF_AMR_KEEP
        endif
      endblock
    enddo

  endsubroutine Estimate_RefinementIndicator2D_t