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.
The estimate runs in two phases. The first is element-local and parallel: it forms each element's modal spectrum, its raw smoothness ratio S_e, and its gate energy g = sum_v w_v E_tot,v. The second applies the amplitude gate, the log10 and the thresholds, and is deferred because the gate needs a field-wide energy scale (a reduction, and under MPI a collective) that no element-local pass can supply. Between the phases indicator(:) transiently holds the raw ratio rather than its logarithm.
comm, when present, is the MPI communicator over which the automatic energy scale is maximized, so that the flags do not depend on the domain decomposition. gate, when present, replaces the weighted-sum gate energy with a caller-computed per-element value - for example an exactly integrated nodal entropy - and is used both for the scale and for the gate comparison.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(RefinementIndicator3D_t), | intent(inout) | :: | this | |||
| class(Scalar3D), | intent(in) | :: | solution | |||
| integer, | intent(in) | :: | ivar | |||
| integer, | intent(in), | optional | :: | comm | ||
| real(kind=prec), | intent(in), | optional | :: | gate(:) |
subroutine Estimate_RefinementIndicator3D_t(this,solution,ivar,comm,gate)
!! 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.
!!
!! The estimate runs in two phases. The first is element-local and parallel: it forms each
!! element's modal spectrum, its raw smoothness ratio S_e, and its gate energy
!! g = sum_v w_v E_tot,v. The second applies the amplitude gate, the log10 and the thresholds,
!! and is deferred because the gate needs a field-wide energy scale (a reduction, and under
!! MPI a collective) that no element-local pass can supply. Between the phases indicator(:)
!! transiently holds the raw ratio rather than its logarithm.
!!
!! comm, when present, is the MPI communicator over which the automatic energy scale is
!! maximized, so that the flags do not depend on the domain decomposition. gate, when present,
!! replaces the weighted-sum gate energy with a caller-computed per-element value - for
!! example an exactly integrated nodal entropy - and is used both for the scale and for the
!! gate comparison.
implicit none
class(RefinementIndicator3D_t),intent(inout) :: this
class(Scalar3D),intent(in) :: solution
integer,intent(in) :: ivar
integer,intent(in),optional :: comm
real(prec),intent(in),optional :: gate(:)
! Local
integer :: iel
real(prec) :: energyScale
call CheckEstimateArguments(this,solution,ivar,gate)
call ResolveEnergyWeights(this,solution%nVar,ivar)
do concurrent(iel=1:this%nElem)
block
integer :: j,k,p,q,rr,ii,v,v0,v1,Np
real(prec) :: tmp(1:this%N+1,1:this%N+1,1:this%N+1)
real(prec) :: uhat(1:this%N+1,1:this%N+1,1:this%N+1)
real(prec) :: acc,etot,eclip1,eclip2,r1,r2,se,semax,g,w
logical :: needSe,needG
! Absolute (machine-epsilon) guard on the energy RATIOS below. This is not the amplitude
! gate: it only keeps 0/0 out of r1 and r2. The amplitude gate is relative and is applied
! in FinalizeIndicator, once the field-wide energy scale is known.
real(prec),parameter :: energyFloor = epsilon(1.0_prec)
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
g = 0.0_prec
do v = 1,solution%nVar
! A variable is transformed if it drives the indicator, or if it carries a non-zero gate
! weight (its energy is then needed for the gate). With the default weights these
! coincide, so the work is exactly what it was before the gate existed.
w = this%energyWeight(v)
needSe = (v >= v0 .and. v <= v1)
needG = (w /= 0.0_prec)
if(.not.(needSe .or. needG)) cycle
! Pass 1 (xi direction): tmp(p,j,k) = sum_i Pmodal(i,p) * u(i,j,k)
do k = 1,Np
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,k,iel,v)
enddo
tmp(p,j,k) = acc
enddo
enddo
enddo
! Pass 2 (eta direction): uhat(p,q,k) = sum_j Pmodal(j,q) * tmp(p,j,k)
do k = 1,Np
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,k)
enddo
uhat(p,q,k) = acc
enddo
enddo
enddo
! Pass 3 (zeta direction): tmp(p,q,r) = sum_k Pmodal(k,r) * uhat(p,q,k), so tmp holds
! the full 3-D modal coefficients uhat(p,q,r) after this pass. The two buffers alternate
! (tmp -> uhat -> tmp) so only two element-sized temporaries are needed, mirroring the
! shared-memory budget of the GPU kernel.
do rr = 1,Np
do q = 1,Np
do p = 1,Np
acc = 0.0_prec
do ii = 1,Np
acc = acc+this%Pmodal(ii,rr)*uhat(p,q,ii)
enddo
tmp(p,q,rr) = acc
enddo
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 rr = 1,Np
do q = 1,Np
do p = 1,Np
etot = etot+tmp(p,q,rr)*tmp(p,q,rr)
if(p <= Np-1 .and. q <= Np-1 .and. rr <= Np-1) then
eclip1 = eclip1+tmp(p,q,rr)*tmp(p,q,rr)
endif
if(p <= Np-2 .and. q <= Np-2 .and. rr <= Np-2) then
eclip2 = eclip2+tmp(p,q,rr)*tmp(p,q,rr)
endif
enddo
enddo
enddo
! Gate energy: g = sum_v w_v E_tot,v, the discrete (quadratic) entropy integral over
! the element when the weights come from an entropy Hessian.
if(needG) g = g+w*etot
if(needSe) then
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)
endif
enddo
! Raw ratio; FinalizeIndicator gates it and takes the log10 in place.
this%indicator(iel) = semax
this%gate(iel) = g
endblock
enddo
if(present(gate)) then
do iel = 1,this%nElem
this%gate(iel) = gate(iel)
enddo
endif
call ResolveEnergyScale(this,energyScale,comm)
call FinalizeIndicator(this,energyScale)
endsubroutine Estimate_RefinementIndicator3D_t