Block based adaptive mesh refinement

[kidder]

Goal: Anisotropic h-p adaptive mesh refinement in D dimensions.

Description of block-based domain

Glossary

• anisotropic refinement: each dimension can refine independently.
• block: the root element of an adaptive-mesh hierarchy.
• direction: the lower or upper side along a coordinate axis.
• domain: a set of non-overlapping blocks.
• element: a D-dimensional mapped reference cell.
• embedding map: a mapping from a d-dimensional coordinate frame to a D-dimensional coordinate frame with d <= D.
• extents: an index specifying the number of grid points along each logical coordinate axis of a manifold.
• grid coordinate frame: a coordinate frame in which the grid is constructed.
• h-refinement: refinement by splitting and joining elements.
• index: a set of integers
• interface: a (D-1)-dimensional mapped reference cell corresponding to either the intersection of two faces of neighboring elements, or the face of an element on an external boundary of the domain.
• logical coordinate frame: the coordinate frame associated with the reference cell.
• manifold: a d-dimensional space embedded in the D-dimensional volume.
• orientation: the relative orientation of a neighboring coordinate system with respect to the local coordinate system.
• p-refinement: refinement by changing the number of grid points in an element.
• physical coordinate frame: a coordinate frame in which the physics of the problem is constructed.
• point: a set of real numbers
• reference cell: a D-dimensional product of the interval [-1,1] along each coordinate axis (i.e. an interval in 1D, a square in 2D, a cube in 3D...).
• regular mesh: a set of discrete grid points on a manifold specified by its extents.
• segment: a subset of the interval [-1,1] along a coordinate axis.

Domain decomposition

• The computational domain is composed of a set of non-overlapping elements.
• The domain is specified from a set of non-overlapping blocks.
• Each face of a block must abut either an external boundary, one other block, or the opposite face of itself (i.e. a periodic boundary in a particular dimension).
• External boundary conditions should not change type within the face of a block.
• Each block is the root element of an adaptive-mesh hierarchy.
• An element is constructed by mapping the reference cell to a grid coordinate frame.
• Each block has a logical coordinate frame and embedding map associated with that of its root element.
• A block is specified by a unique ID, the embedding map of its root element, the directions in which the block abuts an external boundary, and information (ID, direction, and orientation) about its neighboring blocks.
• An element is a subset of a block described by the segments along each coordinate axis of the block.
• An element is specified by a unique ID, the embedding map from the logical coordinate frame of the reference cell to a grid coordinate frame, and information (ID, direction, and orientation) about its neighboring elements.
• The ID of an element is specified by the ID of the block in which it lies, along with the IDs of the segments that specify its location in the block.
• The ID of a segment is specified by the refinement level to which it corresponds, and an index that labels it for the given refinement level. For refinement level K, there are 2^K indices labeling the segments in order from the lower bound of the interval [-1,1] to its upper bound. Each segment on a given refinement level has the same logical coordinate length.

Domain refinement

AMR algorithm

• Refinement can occur whenever a user-specified trigger is triggered.
• User-specified criteria are applied to generate a refinement decision in each dimension of an element.
• A refinement criterion may or may not depend on data from neighboring elements.
• The refinement decision for each direction in an element will be the highest priority of all the specified refinement criteria in that direction for the element.
• Elements can join in a given direction only if they are siblings (i.e. have the same parent segment in that direction), and will have the same refinement level in the transverse directions.

AMR policies

There are a number of potential policies that govern what AMR is allowed to do.

• Allowed type of refinement: h (splitting and joining elements); p (changing resolution in an element); hp (either h or p)
• Isotropic vs Anisotropic refinement (i.e. whether a refinement choice is made for all dimensions, or each dimension individually)
• The priority of refinement choices if multiple refinement criteria are used (e.g. from highest priority to lowest: split, increase resolution, do nothing, decrease resolution, join)
• Allowed difference in refinement level (or resolution) between neighboring elements in the normal (or transverse) direction(s) to the interface between the neighboring elements
• What is done if an invalid decision is made at the bounds of the allowed refinement levels or resolutions

Thoughts on policy choices

• In each transverse direction, h-refinement levels of neighboring elements must be within one refinement level of each other.
• We may enforce a similar restriction for p-refinement, allow user to set? Pass along shock detection information?
• A root segment cannot be joined; instead flag it to decrease resolution.
• An element with extents of two in a given direction cannot decrease resolution; instead flag it to do nothing (or join?).
• A segment at the highest allowed refinement level cannot be split; instead increase resolution (unless near a shock?).
• An element with extents of the allowed maximum in a given direction cannot increase resolution; instead flag it to split; if it cannot be split as it is already at the highest allowed refinement level, abort as either something has gone wrong or you are not resolving the problem.

Development plan

1. Develop domain infrastructure for h-refinement, namely splitting and joining elements, and updating all elements to have correct information about their neighbors.
2. Develop infrastructure for transferring all data to newly created elements, and for (un)registering new (old) elements with other components (e.g. observers, interpolators...)
3. Develop infrastructure for choosing refinement criteria
4. Add remaining infrastructure to enable DG evolution with h-refinement
5. Add additional infrastructure to enable DG-subcell evolution with h-refinement
6. Add p-refinement

Domain infrastructure for h-refinement

This can be worked on with a simple executable that creates an initial domain and then randomly refines elements a specified number of iterations.

• The executable will have two components:
  1. Element: an array component representing the elements
  2. AmrMonitor: a singleton component that gathers diagnostic information in order to check the validity of the refined domain
• The executable will have four phases, executing the listed actions (on the Elements unless otherwise noted):
  1. Initialization: run once to setup the initial domain
     • InitializeDomain: sets up initial domain
     • InitializeAmr: sets up information specific to AMR
  2. CheckDomain: run each iteration to check the validity of the domain
     • SendGlobalAmrDiagnostic: sends global diagnostics to the AmrMonitor via a reduction
     • CheckDomain: run by the AmrMonitor to verify the domain (throws error if Domain is invalid)
  3. EvaluateAmrCriteria: run each iteration to determine how each Element will be refined
     • EvaluateRefinementCriteria: evaluates the refinement criteria to generate an AMR decision on an Element, and sends the decision to the neighbors of the element
     • UpdateAmrDecision: takes the decision of a neighboring Element and determines whether or not the current Element needs to update its AMR decision in order to satisfy AMR policies (e.g. refinement levels need to be within one, not able to join because sibling did not choose to join, ...) If the decision is changed, send the new decision to neighbors.
  4. AdjustDomain: run each iteration to adjust the domain according to the refinement decisions of each element
     • AdjustDomain:
       • in lowest dimension to be split, create the new child Elements and call InitializeChild on each of them, then delete the current element
       • if no dimensions are to be split, determine the best dimension that wants to join with its sibling; then if this Element is the lower sibling of two Elements joining, create the new parent Element; each child Element of the new parent Element calls InitializeParent and then deletes itself
       • if the Element does not want to join or split it updates its neighbor information.

Initialization after refinement

Here are the current tags that need to be mutated for scalar wave:

  // mutable creation tags added by DistributedObject
  using distributed_object_tags =
      tmpl::list<Tags::MetavariablesImpl<Metavariables>,
                 Tags::ArrayIndexImpl<ElementId<volume_dim>>,
                 Tags::GlobalCacheProxy<Metavariables>>;
  // These tags are all initialized when a distributed object is created,
  // so nothing needs to be done for them during refinement

  // mutable creation tags added by evolution::dg::Initialization::Domain
  using domain_tags = tmpl::list<
      ::domain::Tags::InitialExtents<Dim>,
      ::domain::Tags::InitialRefinementLevels<Dim>,
      evolution::dg::Tags::Quadrature, evolution::dg::Tags::NeighborMesh<Dim>,
      ::domain::Tags::Mesh<Dim>, ::domain::Tags::Element<Dim>,
      ::domain::Tags::ElementMap<Dim, Frame::Grid>,
      ::domain::CoordinateMaps::Tags::CoordinateMap<Dim, Frame::Grid,
                                                    Frame::Inertial>>;
  // only Mesh, NeighborMesh and Element need to be updated

  // mutable creation tags added by Initialization::TimeStepping
  using TimeStepperType =
      tmpl::conditional_t<UsingLts, LtsTimeStepper, TimeStepper>;
  using common_time_stepping_tags =
      tmpl::list<::Tags::Next<::Tags::TimeStepId>, ::Tags::TimeStep,
                 ::Tags::Next<::Tags::TimeStep>, ::Tags::TimeStepId,
                 ::Tags::AdaptiveSteppingDiagnostics, ::Tags::Time,
                 Tags::InitialTimeDelta, Tags::InitialSlabSize<UsingLts>,
                 ::Tags::TimeStepper<TimeStepperType>>;
  using conditional_lts_tags = tmpl::list<::Tags::StepChoosers>;
  using step_chooser_tags =
      StepChoosers::step_chooser_simple_tags<Metavariables, UsingLts>;
  using error_control_tags =
      tmpl::list<::Tags::StepperError<EvolvedVariableTag>,
                 ::Tags::PreviousStepperError<EvolvedVariableTag>,
                 ::Tags::StepperErrorUpdated>;
  static_assert(std::is_same_v<step_chooser_tags, tmpl::list<>> or
                    std::is_same_v<step_chooser_tags, error_control_tags>,
                "step_chooser_tags contains a tag that is not handled for amr");
  // error_contorl_tags need to be projected

  // mutable creation tags added by Initialization::TimeStepperHistory
  using variables_tag = typename Metavariables::system::variables_tag;
  using dt_variables_tag = db::add_tag_prefix<::Tags::dt, variables_tag>;
  using time_stepper_history_tags =
      tmpl::list<dt_variables_tag,
                 ::Tags::HistoryEvolvedVariables<variables_tag>>;
  // time_stepper_history_tags need to be projected

  // mutable creation tags added by
  // Initialization::Actions::NonconservativeSystem
  using nonconservative_sytem_tags = tmpl::list<variables_tag>;
  // variables_tag needs to be projected

  // mutable creation tags added by
  // evolution::Initialization::Actions::SetVariables
  using set_variables_tags = tmpl::list<::Tags::Time>;
  // this tag was already taken care of...

  // mutable creation tags added by ScalarWave::Actions::InitializeConstraints
  using constraint_tags = tmpl::list<ScalarWave::Tags::ConstraintGamma2>;
  // constraint_tag needs to be projected

  // mutable creation tags added by interp::Actions::ElementInitInterpPoints
  using interp_tags = tmpl::list<InterpPointInfoTag>;
  // ????

  // mutable creation tags added by ::evolution::dg::Initialization::Mortars
  using Dim = Metavariables::volume_dim;
  using System = typename Metavariables::system;
  using mortar_tags = tmpl::list<
      ::domain::Tags::InitialExtents<Dim>, evolution::dg::Tags::Quadrature,
      Tags::MortarData<Dim>, Tags::MortarMesh<Dim>, Tags::MortarSize<Dim>,
      Tags::MortarNextTemporalId<Dim>,
      evolution::dg::Tags::NormalCovectorAndMagnitude<Dim>,
      Tags::MortarDataHistory<
          Dim, typename db::add_tag_prefix<
                   ::Tags::dt, typename System::variables_tag>::type>>;
  // mortar_tags either need to be updated, projected, or cleared

  // mutable creation tags added by
  // evolution::Actions::InitializeRunEventsAndDenseTriggers
  using event_and_dense_trigger_tags =
      tmpl::list<evolution::Tags::EventsAndDenseTriggers,
                 Tags::PreviousTriggerTime>;
  // should be okay as is

  // mutable creation tags add by SelfStart::Actions::Initialize
  using self_start_tags =
      tmpl::list<SelfStart::Tags::InitialValue<Tags::TimeStep>,
                 SelfStart::Tags::InitialValue<Tags::Next<Tags::TimeStep>>,
                 SelfStart::Tags::InitialValue<variables_tag>>;
  // should be okay as is

[nilsvu]

Q: When an element is split, do the two new elements retain the p refinement of their parent (effectively doubling the number of grid points), or decrease their p refinement?

[kidder]

That could be a policy. Initially, the new elements will retain their p-refinement.