SELF_RefinementIndicator_3D_t Module

  1. SELF_RefinementIndicator_3D_t.f90
  2. SELF_RefinementIndicator_3D_t

Modal (Legendre spectral-decay) refinement indicator for 3-D spectral element solutions.

Purpose

Provides the trigger mechanism for adaptive mesh refinement (AMR): a per-element scalar that measures how much of a solution field's energy sits in its highest polynomial modes. A well-resolved (smooth) field has energy that decays rapidly with mode number, so almost all of its energy lives in the low modes; an under-resolved field (a steep front, a discontinuity, or simply too-coarse an element) leaves significant energy in the top modes. Comparing the top-mode energy fraction against thresholds flags each element for refinement, coarsening, or no change.

Method (Legendre modal-energy / spectral-decay indicator)

Following Persson & Peraire (2006), "Sub-cell shock capturing for discontinuous Galerkin methods", AIAA 2006-112, and the robustified variant of Hennemann, Rueda-Ramirez, Hindenlang & Gassner (2021), J. Comput. Phys. 426, 109935 (the shock indicator used in Trixi.jl), the nodal solution on an element is expanded in a tensor-product Legendre modal basis

u(xi,eta,zeta) = sum_{p=0}^{N} sum_{q=0}^{N} sum_{r=0}^{N} uhat(p,q,r) Ltilde_p(xi) Ltilde_q(eta) Ltilde_r(zeta),

where Ltilde_p are the L2-normalized Legendre polynomials on [-1,1] (int_{-1}^{1} Ltilde_p Ltilde_q dxi = delta_pq). With that normalization the total modal energy equals the exact L2 energy of the field on the reference element,

E_tot = sum_{p,q,r} uhat(p,q,r)^2 = ||u||_{L2([-1,1]^3)}^2 .

Two "clipped" energies are formed by dropping the highest one and two modes in each direction,

E_clip1 = sum_{p,q,r <= N-1} uhat(p,q,r)^2 , E_clip2 = sum_{p,q,r <= N-2} uhat(p,q,r)^2 ,

and the smoothness ratio is the larger of the top-shell and next-shell energy fractions,

S_e = max( (E_tot - E_clip1)/E_tot , (E_clip1 - E_clip2)/E_clip1 ) .

The second term (present only for N >= 2) guards against the odd/even parity dropouts that a single-mode measure can suffer for symmetric data - this is the Hennemann-Gassner refinement of the original single-mode Persson-Peraire estimate. The indicator returned per element is the base-10 logarithm sigma_e = log10(S_e); a near-machine-zero floor keeps it finite for a perfectly resolved (or identically zero) field.

Amplitude gate (relative energy floor)

S_e is a ratio and therefore scale-free: it says nothing about how much energy an element carries, only how that energy is distributed across modes. An element holding a negligible share of the field's energy - low-amplitude grid-scale residue in the wake of a passing wave, say - still shows a flat modal spectrum and so reads as under-resolved forever. Guarding the ratio with an absolute floor at machine epsilon does not help: E_tot carries the SQUARE of the field amplitude, so a wake at 1e-5 of the peak amplitude sits ~1e10 above epsilon and still takes the ratio branch. Such elements are never flagged for coarsening and the mesh behind a front is never released.

The fix is a second, RELATIVE floor. A per-element gate energy

g_e = sum_v w_v E_tot,e,v

is formed from the modal energies of the driving variables with non-negative weights w_v. Since E_tot,e,v is the exact L2 energy of variable v on the reference element, g_e with w taken from a quadratic entropy function is the discrete entropy integral over the element: a convex function of the state, and therefore an amplitude measure that is meaningful across variables of very different magnitude. Comparing it against a field-wide energy scale,

effective floor = max( epsilon, relativeEnergyFloor * energyScale ) g_e <= effective floor -> S_e := 0 (element is quiescent, hence resolved)

gives an amplitude gate that is independent of the modal shape. Because energy goes as amplitude squared, relativeEnergyFloor is 10**(dB/10) in amplitude terms: the default 1e-12 treats amplitudes below 1e-6 (-120 dB) of the field scale as quiescent.

Raising the floor is NOT free, and the useful range is narrower than it looks. The gate cannot distinguish residue from the low-amplitude flank of a feature that is genuinely under-resolved, so an aggressive floor also suppresses refinement DEPTH: on the ultrasound benchmark a floor of 1e-8 stopped the source pulse from reaching the level cap, which doubled the level-based time step and produced a mesh roughly twice as large by the end of the run. Measure before raising it; see the note on SELF_AMR_DEFAULT_RELFLOOR.

The gate is a single hard cut by default, which reintroduces on the energy axis exactly the thrashing the two sigma thresholds exist to prevent: an element whose energy drifts across that one value flips COARSEN <-> REFINE on successive epochs. Supplying significantEnergyFloor opens a hysteresis band instead,

g_e <= quiescent floor -> COARSEN, whatever the spectrum says quiescent < g_e <= significant -> SELF_AMR_KEEP (hold the mesh) g_e > significant floor -> the spectrum decides, as before

the middle zone reading "too weak to be worth spending levels on, too strong to declare resolved". Inside the band sigma_e still reports the true spectrum; only the flag is held.

energyScale defaults to the largest gate energy over the elements (globally reduced when Estimate is given an MPI communicator), which makes the gate track a decaying front. It can be pinned with SetEnergyScale for callers that need the flag decision to be independent of a floating-point reduction, and the gate itself can be supplied outright through the optional gate argument to Estimate (for example an exactly integrated nodal entropy).

Note that raising coarsenThreshold is NOT a substitute: in a low-amplitude wake sigma_e is near 0, so any threshold high enough to coarsen the wake also stops a genuine front from refining. The two decisions cannot be separated by thresholds alone.

Trigger semantics

With user-supplied thresholds sigma_refine > sigma_coarsen the per-element flag is

sigma_e > sigma_refine -> SELF_AMR_REFINE (+1) top modes carry too much energy sigma_e < sigma_coarsen -> SELF_AMR_COARSEN (-1) field is over-resolved on this element otherwise -> SELF_AMR_KEEP (0)

Larger (closer to zero) sigma_e means a less smooth / less resolved field. Recommended starting values for double precision are sigma_refine ~ -3.0 and sigma_coarsen ~ -8.0; they are problem dependent and are deliberately left as required arguments rather than hidden defaults.

Units and input ranges

  • The indicator is dimensionless (an energy fraction); the driving field may carry any units.
  • Requires N >= 1 (a degree-0 element has no modal spectrum to measure).
  • The nodal->modal transform is the exact inverse of the Legendre Vandermonde built from the interpolant control points, so the indicator is independent of the control-node type (Legendre-Gauss or Legendre-Gauss-Lobatto) and is exact for polynomial data.

This module contains the portable (CPU) implementation using do concurrent over elements. The backend extension modules (SELF_RefinementIndicator_3D) add device storage and a GPU kernel while preserving the mathematics bit-for-bit in structure.


Uses

  • module~~self_refinementindicator_3d_t~~UsesGraph module~self_refinementindicator_3d_t SELF_RefinementIndicator_3D_t module~self_lagrange~2 SELF_Lagrange module~self_refinementindicator_3d_t->module~self_lagrange~2 mpi mpi module~self_refinementindicator_3d_t->mpi module~self_constants SELF_Constants module~self_refinementindicator_3d_t->module~self_constants iso_c_binding iso_c_binding module~self_refinementindicator_3d_t->iso_c_binding module~self_scalar_3d~2 SELF_Scalar_3D module~self_refinementindicator_3d_t->module~self_scalar_3d~2 module~self_lagrange~2->module~self_constants module~self_lagrange~2->iso_c_binding iso_fortran_env iso_fortran_env module~self_lagrange~2->iso_fortran_env module~self_lagrange_t SELF_Lagrange_t module~self_lagrange~2->module~self_lagrange_t module~self_constants->iso_c_binding module~self_constants->iso_fortran_env module~self_scalar_3d~2->module~self_constants module~self_scalar_3d~2->iso_c_binding module~self_scalar_3d_t SELF_Scalar_3D_t module~self_scalar_3d~2->module~self_scalar_3d_t module~self_scalar_3d_t->module~self_lagrange~2 module~self_scalar_3d_t->module~self_constants module~self_scalar_3d_t->iso_c_binding module~self_datapool SELF_DataPool module~self_scalar_3d_t->module~self_datapool module~self_hdf5 SELF_HDF5 module~self_scalar_3d_t->module~self_hdf5 HDF5 HDF5 module~self_scalar_3d_t->HDF5 module~self_data SELF_Data module~self_scalar_3d_t->module~self_data FEQParse FEQParse module~self_scalar_3d_t->FEQParse module~self_metadata SELF_Metadata module~self_scalar_3d_t->module~self_metadata module~self_lagrange_t->module~self_constants module~self_lagrange_t->iso_c_binding module~self_lagrange_t->iso_fortran_env module~self_lagrange_t->module~self_hdf5 module~self_lagrange_t->HDF5 module~self_quadrature SELF_Quadrature module~self_lagrange_t->module~self_quadrature module~self_supportroutines SELF_SupportRoutines module~self_lagrange_t->module~self_supportroutines module~self_datapool->module~self_constants module~self_hdf5->mpi module~self_hdf5->module~self_constants module~self_hdf5->iso_fortran_env module~self_hdf5->HDF5 module~self_data->module~self_lagrange~2 module~self_data->module~self_constants module~self_data->iso_c_binding module~self_data->module~self_hdf5 module~self_data->HDF5 module~self_data->FEQParse module~self_data->module~self_metadata module~self_metadata->module~self_hdf5 module~self_metadata->HDF5 module~self_quadrature->module~self_constants module~self_quadrature->iso_fortran_env module~self_supportroutines->module~self_constants module~self_supportroutines->iso_fortran_env

Used by

  • module~~self_refinementindicator_3d_t~~UsedByGraph module~self_refinementindicator_3d_t SELF_RefinementIndicator_3D_t module~self_refinementindicator_3d SELF_RefinementIndicator_3D module~self_refinementindicator_3d->module~self_refinementindicator_3d_t module~self_refinementindicator_3d~2 SELF_RefinementIndicator_3D module~self_refinementindicator_3d~2->module~self_refinementindicator_3d_t module~self_amrcontroller_3d SELF_AMRController_3D module~self_amrcontroller_3d->module~self_refinementindicator_3d

Contents

Variables

Derived Types

Functions

Subroutines


Variables

TypeVisibilityAttributesNameInitial
integer, public, parameter:: SELF_AMR_ALLVARS =0
integer, public, parameter:: SELF_AMR_COARSEN =-1
real(kind=prec), public, parameter:: SELF_AMR_DEFAULT_RELFLOOR =1.0e-12_prec
integer, public, parameter:: SELF_AMR_KEEP =0
integer, public, parameter:: SELF_AMR_REFINE =1

Derived Types

type, public :: RefinementIndicator3D_t

Components

TypeVisibilityAttributesNameInitial
integer, public :: N =0

Polynomial degree of the interpolant the indicator is built for.
real(kind=prec), public, pointer, contiguous, dimension(:,:):: Pmodal=> null()

Nodal-to-modal transform. Pmodal(ii,p) is the (p,ii) entry of the inverse Legendre Vandermonde in the L2-normalized basis, so the 1-D modal coefficients are uhat(p) = sum_ii Pmodal(ii,p) * u(ii) (first index summed, SELF matrix convention).
real(kind=prec), public :: coarsenThreshold =0.0_prec

Elements with sigma_e below this value are flagged SELF_AMR_COARSEN.
real(kind=prec), public :: energyScale =0.0_prec

Squared field scale the relative floor is measured against, in the units of the gate energy. Meaningful only when energyScaleIsSet is true; otherwise the scale is recomputed from the current field on every Estimate.
logical, public :: energyScaleIsSet =.false.

Whether energyScale was pinned by SetEnergyScale (true) or is computed automatically as the largest gate energy over the elements (false, the default).
real(kind=prec), public, pointer, contiguous, dimension(:):: energyWeight=> null()

Non-negative weights w_v of the gate energy g_e = sum_v w_v E_tot,e,v, indexed by solution variable. Resolved lazily, because the variable count is a property of the solution field and is not known at Init.
logical, public :: energyWeightsSet =.false.

Whether energyWeight was supplied by SetEnergyWeights (true) or is regenerated from the driving-variable index on every Estimate (false, the default).
integer, public, pointer, contiguous, dimension(:):: flag=> null()

Per-element refinement flag: SELF_AMR_REFINE / SELF_AMR_KEEP / SELF_AMR_COARSEN.
real(kind=prec), public, pointer, contiguous, dimension(:):: gate=> null()

Per-element gate energy g_e from the most recent Estimate. Retained as a diagnostic: it is the quantity the amplitude gate actually compared against the effective floor.
real(kind=prec), public, pointer, contiguous, dimension(:):: indicator=> null()

Per-element indicator value sigma_e = log10(S_e).
integer, public :: nElem =0

Number of (rank-local) elements the indicator arrays are sized for.
integer, public :: nVarWeights =0

Allocated length of energyWeight (the solution variable count it was resolved for).
real(kind=prec), public :: refineThreshold =0.0_prec

Elements with sigma_e above this value are flagged SELF_AMR_REFINE.
real(kind=prec), public :: relativeEnergyFloor =SELF_AMR_DEFAULT_RELFLOOR

Elements whose gate energy is at or below relativeEnergyFloorenergyScale are treated as quiescent (hence perfectly resolved) regardless of their modal shape. Energy goes as amplitude squared, so this is 10*(dB/10) in amplitude terms: 1e-8 gates amplitudes below 1e-4 (-80 dB) of the field scale. Set to 0 to recover the pure absolute (machine-epsilon) floor.
real(kind=prec), public :: significantEnergyFloor =SELF_AMR_DEFAULT_RELFLOOR

Upper edge of the hysteresis band on the energy axis. Elements between relativeEnergyFloor and this fraction of the energy scale are flagged SELF_AMR_KEEP whatever their spectrum: too weak to justify spending levels on, too strong to declare resolved. Equal to relativeEnergyFloor by default, which collapses the band to a single hard cut. See SetRelativeEnergyFloor.

Type-Bound Procedures

procedure, public :: ClearEnergyScale => ClearEnergyScale_RefinementIndicator3D_t
procedure, public :: CountFlagged => CountFlagged_RefinementIndicator3D_t
procedure, public :: Estimate => Estimate_RefinementIndicator3D_t
procedure, public :: Free => Free_RefinementIndicator3D_t
procedure, public :: Init => Init_RefinementIndicator3D_t
procedure, public :: SetEnergyScale => SetEnergyScale_RefinementIndicator3D_t
procedure, public :: SetEnergyWeights => SetEnergyWeights_RefinementIndicator3D_t
procedure, public :: SetRelativeEnergyFloor => SetRelativeEnergyFloor_RefinementIndicator3D_t
procedure, public :: SetThresholds => SetThresholds_RefinementIndicator3D_t
procedure, public :: UpdateDevice => UpdateDevice_RefinementIndicator3D_t
procedure, public :: UpdateHost => UpdateHost_RefinementIndicator3D_t

Functions

public function CountFlagged_RefinementIndicator3D_t(this, flagValue) result(n)

Count the rank-local elements currently carrying the requested flag value (SELF_AMR_REFINE / SELF_AMR_KEEP / SELF_AMR_COARSEN). Provided as a convenience for drivers and tests; a global count across MPI ranks is the caller's responsibility.

Arguments

TypeIntentOptionalAttributesName
class(RefinementIndicator3D_t), intent(in) :: this
integer, intent(in) :: flagValue

Return Value integer

public pure function NormalizedLegendre(p, x) result(Lp)

L2-normalized Legendre polynomial Ltilde_p(x) = L_p(x)*sqrt((2p+1)/2) on [-1,1], evaluated with the standard three-term recurrence. The normalization gives int_{-1}^{1} Ltilde_p Ltilde_q dx = delta_pq.

Arguments

TypeIntentOptionalAttributesName
integer, intent(in) :: p
real(kind=real64), intent(in) :: x

Return Value real(kind=real64)


Subroutines

public subroutine BuildModalTransform(controlPoints, N, Pmodal)

Build the nodal->modal transform Pmodal for a 1-D degree-N interpolant whose nodes are controlPoints(1:N+1). The transform is the exact inverse of the L2-normalized Legendre Vandermonde V(i,p) = Ltilde_{p-1}(x_i); Pmodal(ii,p) stores V^{-1}(p,ii) so that uhat(p) = sum_ii Pmodal(ii,p) * u(ii). The Vandermonde inverse is formed in double precision for conditioning and cast back to the working precision.

Arguments

TypeIntentOptionalAttributesName
real(kind=prec), intent(in) :: controlPoints(1:N+1)
integer, intent(in) :: N
real(kind=prec), intent(out) :: Pmodal(1:N+1,1:N+1)

public subroutine CheckEstimateArguments(this, solution, ivar, gate)

Argument validation shared by the portable and backend Estimate implementations.

Arguments

TypeIntentOptionalAttributesName
class(RefinementIndicator3D_t), intent(in) :: this
class(Scalar3D), intent(in) :: solution
integer, intent(in) :: ivar
real(kind=prec), intent(in), optional :: gate(:)

public subroutine ClearEnergyScale_RefinementIndicator3D_t(this)

Return to the automatic energy scale (the largest gate energy over the elements).

Arguments

TypeIntentOptionalAttributesName
class(RefinementIndicator3D_t), intent(inout) :: this

public 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.

Read more…

Arguments

TypeIntentOptionalAttributesName
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(:)

public subroutine FinalizeIndicator(this, energyScale)

Second phase of an estimate, shared by every backend so that the flags are identical whichever computed the spectra: apply the amplitude gate, take the log10 and set the refine/keep/coarsen flags.

Read more…

Arguments

TypeIntentOptionalAttributesName
class(RefinementIndicator3D_t), intent(inout) :: this
real(kind=prec), intent(in) :: energyScale

public subroutine Free_RefinementIndicator3D_t(this)

Arguments

TypeIntentOptionalAttributesName
class(RefinementIndicator3D_t), intent(inout) :: this

public subroutine Init_RefinementIndicator3D_t(this, interp, nElem, refineThreshold, coarsenThreshold)

Allocate the indicator for an interpolant of degree interp%N and nElem elements and precompute the nodal->modal transform matrix from the interpolant control points.

Read more…

Arguments

TypeIntentOptionalAttributesName
class(RefinementIndicator3D_t), intent(out) :: this
type(Lagrange), intent(in), target:: interp
integer, intent(in) :: nElem
real(kind=prec), intent(in) :: refineThreshold
real(kind=prec), intent(in) :: coarsenThreshold

public subroutine InvertMatrix(A, Ainv, n)

Invert the n-by-n matrix A by Gauss-Jordan elimination with partial pivoting, in double precision. Used once at initialization for the small (N+1) Legendre Vandermonde; A is a well-conditioned Vandermonde in an orthonormal basis, so a direct solve is appropriate.

Arguments

TypeIntentOptionalAttributesName
real(kind=real64), intent(in) :: A(1:n,1:n)
real(kind=real64), intent(out) :: Ainv(1:n,1:n)
integer, intent(in) :: n

public subroutine ResolveEnergyScale(this, energyScale, comm)

Determine the energy scale the relative floor is measured against: the pinned value when SetEnergyScale was used, otherwise the largest gate energy over the (rank-local) elements. When comm is present the maximum is taken over the communicator, so every rank applies the same floor and the flags - hence the adapted mesh - do not depend on the decomposition. One small collective per indicator evaluation, i.e. per adaptation epoch, never inside the time-stepping loop.

Arguments

TypeIntentOptionalAttributesName
class(RefinementIndicator3D_t), intent(in) :: this
real(kind=prec), intent(out) :: energyScale
integer, intent(in), optional :: comm

public subroutine ResolveEnergyWeights(this, nVar, ivar)

Ensure energyWeight is associated and sized to nVar. Weights supplied by SetEnergyWeights must match nVar; otherwise the default weights are regenerated: unit weight on the driving variable, or on every variable when ivar is SELF_AMR_ALLVARS. Resolved here rather than at Init because the variable count belongs to the solution field, not to the indicator.

Arguments

TypeIntentOptionalAttributesName
class(RefinementIndicator3D_t), intent(inout) :: this
integer, intent(in) :: nVar
integer, intent(in) :: ivar

public subroutine SetEnergyScale_RefinementIndicator3D_t(this, energyScale)

Pin the energy scale that normalizes the relative floor, in the units of the gate energy (squared field amplitude times the reference-element volume). Two reasons to use this rather than the automatic maximum over elements:

Read more…

Arguments

TypeIntentOptionalAttributesName
class(RefinementIndicator3D_t), intent(inout) :: this
real(kind=prec), intent(in) :: energyScale

public subroutine SetEnergyWeights_RefinementIndicator3D_t(this, w)

Set the per-variable weights w_v >= 0 of the gate energy g_e = sum_v w_v E_tot,e,v, where E_tot,e,v is the exact L2 energy of variable v on the reference element.

Read more…

Arguments

TypeIntentOptionalAttributesName
class(RefinementIndicator3D_t), intent(inout) :: this
real(kind=prec), intent(in) :: w(:)

public subroutine SetRelativeEnergyFloor_RefinementIndicator3D_t(this, relativeEnergyFloor, significantEnergyFloor)

Set the relative energy floor of the amplitude gate. An element whose gate energy g satisfies

Read more…

Arguments

TypeIntentOptionalAttributesName
class(RefinementIndicator3D_t), intent(inout) :: this
real(kind=prec), intent(in) :: relativeEnergyFloor
real(kind=prec), intent(in), optional :: significantEnergyFloor

public subroutine SetThresholds_RefinementIndicator3D_t(this, refineThreshold, coarsenThreshold)

Update the refine/coarsen thresholds without rebuilding the transform.

Arguments

TypeIntentOptionalAttributesName
class(RefinementIndicator3D_t), intent(inout) :: this
real(kind=prec), intent(in) :: refineThreshold
real(kind=prec), intent(in) :: coarsenThreshold

public subroutine UpdateDevice_RefinementIndicator3D_t(this)

Arguments

TypeIntentOptionalAttributesName
class(RefinementIndicator3D_t), intent(inout) :: this

public subroutine UpdateHost_RefinementIndicator3D_t(this)

Arguments

TypeIntentOptionalAttributesName
class(RefinementIndicator3D_t), intent(inout) :: this