\section{INTRODUCTION}

Many problems in astrophysics and cosmology will press against the boundaries of supercomputer software and hardware development for some time to come. 
Among the most challenging are time-dependent systems that should be treated in three position space dimensions, plus up to three momentum space dimensions for those that involve radiation transport. 
Several types of physics---and their couplings---must be addressed. 
Taken together, physics such as self-gravity, turbulent cascades, reactions that may or may not be in equilibrium, and the operation of both microscopic and macroscopic processes often imply the simultaneous relevance of multiple scales in space and time, typically requiring some kind of spatial adaptivity and multiple solvers (at least some of which may be time-implicit). 
Solvers deployed with good resolution in as much of phase space as possible fill the memory and churn through the cycles of the largest supercomputers available. 
That the relevant theory is described by time-dependent partial (and sometimes integro-)differential equations implies the need for synchronous evolution and communication between different regions of position space (and sometimes momentum space). 
These latter aspects seem ever more difficult to address due to the fact that increases in computing capability seem only to come with such additional burdens as distributed memory and distributed (and, more recently, heterogeneous) processing capacity. 
From the perspective of working astrophysicists and cosmologists, it can seem that the physics itself recedes ever further away as their efforts are channeled towards developing and tailoring codes to the specific features of these increasingly complex high-performance machines. 
In this environment, the availability of well-designed codes with broadly applicable physics capabilities is increasingly valuable to researchers.

One such computationally demanding problem is the elucidation of the explosion mechanism of core-collapse supernovae \citep{Mezzacappa2005,Woosley2005,Kotake2006,JANKA2007,Janka2012}. 
A star with mass greater than $\sim 8\, M_\odot$ develops a degenerate core during its final burning stage. 
Once the core becomes sufficiently large---roughly the Chandrasekhar mass---it becomes unstable and undergoes catastrophic collapse on a nearly free-fall time scale. 
Collapse of an inner subsonic portion of the core halts around nuclear density (where nucleons begin to overlap), leading to the development of a shock wave at the interface between this inner core and the supersonically infalling outer portion of the core. 
The shock eventually will disrupt the entire star and give rise to the luminous supernova, but it stagnates shortly after its formation due to the endothermic reduction of heavy nuclei into their constituent nucleons and enervating neutrino losses. 
The mechanism of shock revival---that is, the explosion mechanism---remains to be fully elucidated by numerical simulations, but is expected to involve some combination of heating by intense neutrino fluxes streaming from the nascent neutron star, fluid instabilities, rotation, and magnetic fields. 
In fact, the relative contributions of these phenomena may vary from event to event: pre-supernova stars differ in properties such as mass and rotation, leading to explosions that may or may not be jet-like with an associated gamma ray burst, and either a black hole or a neutron star as a compact remnant.

Even this brief description conveys an initial sense of the multiphysics and multiscale nature of core-collapse supernova explosion mechanism simulations. 
The problem manifestly requires self-gravity, which must be general relativistic to treat the more rare and extreme case of black hole formation, and ideally would be at least approximately relativistic even when the compact remnant is a neutron star. 
The treatment of hydrodynamics---and ideally magnetohydrodynamics, especially in connection with black hole formation---must be able to handle shocks. At high density the equation of state must describe neutron-rich nuclear matter at finite temperature, and at low density it is desirable to track nuclear composition with a reaction network spanning a wide range of nuclear species. 
Neutrino transport must span diffusive, decoupling, and free-streaming regimes, and include several species and their interactions with each other and with various fluid constituents. 
Gravitational collapse, a steepening density cliff at the surface of the nascent neutron star, and regions of turbulence strongly recommend some form of spatial adaptivity.
The stiff equations of neutrino transport and a nuclear network normally require time-implicit evolution. 
At the present time computational limitations still require dimensional reduction of phase space; simple estimates suggest that at least exascale resources will be required for the full transport problem. 
There is fairly wide agreement that retention of at least energy dependence in full neutrino transport is important, and that all three position space dimensions should be included, but finished simulations of this type have yet to be published by any group \citep[though some are in progress; see for instance][]{Bruenn2009}.

GenASiS\ ({\em Gen}eral {\em A}strophysical {\em Si}mulation {\em S}ystem) is a new code being developed, at least initially and primarily, for the simulation of core-collapse supernovae on the world's leading capability supercomputers. 
`General' denotes the capacity of the code to include and refer to multiple algorithms, solvers, and physics and numerics choices with the same abstracted names and/or interfaces. 
In GenASiS\ this is accomplished with features of Fortran 2003 that support the object-oriented programming paradigm \citep[e.g.][]{Reid2007,Adams2008}.
`Astrophysical' roughly suggests---over-broadly, at least initially---the types of systems at which the code is aimed, and the kinds of physics and solvers it makes available. 
`Simulation System' indicates that the code is not a single program, but a collection of modules, structured as classes, that can be invoked by a suitable driver program set up to characterize and initialize a particular problem.

Given that several codes used for astrophysical simulations in general and core-collapse supernova simulations in particular have been described in the literature, an ambitious, new, from-the-ground-up undertaking like GenASiS\ should offer something new in terms of approach, capability, and/or availability. 

The potential generality, extensibility, and maintainability afforded by an object-oriented approach is one such distinguishing feature.  
This programming paradigm involves several interrelated principles.
{\em Abstraction} identifies the major concepts required by a program, without specifying the details of implementation.
{\em Encapsulation} bundles together and controls access to the data (or {\em members}) and actions (or {\em methods}) associated with a particular concept into a self-contained unit---a {\em class}.
(The class itself is just a kind of template or glorified type specification; an {\em instance} of the class, basically a variable declared with that type specification, is called an {\em object}.
In terms of an analogy with intrinsic data types, a class is to an object as the intrinsic data type {\tt real} is to a variable declared as {\tt real}.)
Done well, abstraction and encapsulation lead to {\em decoupling}---the separation of code into building blocks that are as independent and reusable as possible. 
Reusability is further enhanced by {\em polymorphism}, which enables multiple versions or implementations of the same basic concept to be referred to and used interchangeably.
Closely related is {\em inheritance,} which allows new ({\em child}) classes to be formed from prior ({\em parent}) classes: existing members and methods are retained or modified, and new members and methods can be added.

Several challenges facing the development and ongoing relevance of complicated simulation codes can be ameliorated by these principles of object-oriented design, and can do so without sacrificing performance.
For a problem as rich as core-collapse supernovae, simulating all the desired physics in full dimensionality right off the bat is not feasible; a series of approximations of increasing sophistication must be deployed over time. (Consider for example Newtonian vs. relativistic gravity, mean heavy nucleus composition vs. networks of nuclei, variable Eddington tensor---with any number of closure schemes---vs. Boltzmann transport, and so on.)
Even within a particular physics approximation, there are many ways to formulate and solve it. (Imagine different coordinate systems, linear and nonlinear solvers, hydrodynamics reconstruction schemes and Riemann solvers, and so forth.) 
Multiple hardware platforms (workstations vs. basic clusters vs. increasingly sophisticated supercomputers with heterogeneous, multicore, accelerator-enhanced nodes), and different software libraries addressing the same needs (e.g. I/O, FFT, solvers, etc.), are available or may come into existence.
Some algorithms may be greatly facilitated with the availability of multiple copies of the same basic entities with their own self-contained storage (think for instance of the multiple levels of both meshing and sets of related variables needed by an adaptive mesh refinement scheme---a major motivation for us, as discussed below).  
All of these challenges are ripe for attack with abstraction, encapsulation, decoupling, polymorphism, and inheritance.
Indeed, it is the object-oriented approach that best allows for the flexibility connoted by the `General' in GenASiS---the capacity of the code to include and refer to multiple algorithms, solvers, and physics and numerics choices with the same abstracted names and/or interfaces.
Moreover, higher-level organization and flexibility afforded by sophisticated data structures need not sacrifice performance.
Lower-level kernels in which the bulk of computational time is spent can be kept isolated and simple---and therefore amenable to compiler optimization---by being written in terms of elementary loops over operations on array variables of intrinsic data types.

Beyond the programming paradigm, one fundamental characteristic of a simulation code that underpins almost everything else is the nature---or even existence, when one considers particle methods---of its meshing, as this is the stage upon which the physics plays out. 
Smoothed-particle hydrodynamics has been widely used in the broader astrophysics community as an efficient way to get to three dimensions, but its accuracy has been controversial \citep[e.g.][]{Agertz2006,Price2008,Springel2010}.  
Unstructured meshes are useful for the complex geometries of engineering contexts, but their high overhead is probably not justified for the more simple geometries of most astrophysics problems.
(On the other hand, the relatively new use of adaptive Voronoi tesselations \citep{Springel2010a} is an interesting new approach that merits further consideration.)   
Moving patches may be useful for following compact objects in orbit \citep[e.g.][]{Scheel2006}, but introduce additional complicated source terms associated with non-inertial reference frames, and are not an obvious fit for more centrally-condensed problems like core-collapse supernovae.
When more structured grid-based approaches are considered, one major choice is whether to use adaptive mesh refinement (AMR) in an effort to deploy computational resources only where needed, and if so, what type of AMR should be used. 

Most existing core-collapse supernova codes use strategies other than AMR to handle gravitational collapse. 
One code uses smoothed-particle hydrodynamics \citep{Fryer2006}. Others use an at least relatively high resolution radial mesh, usually in only one \citep{Rampp2002,Thompson2003,Liebendorfer2004,Sumiyoshi2005} or two \citep{Buras2006a,Livne2007,Bruenn2009} position space dimensions, often with Lagrangian coordinates (in spherical symmetry) or a moving radial mesh to follow the infall. 
Use of a radial mesh and associated spherical coordinates imposes severe time step limitations at coordinate singularities unless special measures are taken, such as use of an unstructured mesh to morph to a different coordinate system at low radius \citep{Livne2007}, or an overlap of a radial mesh with a Cartesian mesh at low radius \citep{Scheidegger2008}, or an overlap of two separate radial meshes in a so-called `yin-yang' configuration \citep{Wongwathanarat2010b}. 
While these approaches may give satisfactory gravitational collapse, with a moving radial mesh also reasonably handling the steepening and moving density cliff at the neutron star surface, they cannot address regions of turbulence as well as AMR can. 
Those codes that do use (or intend to use) AMR use the block-structured variety \citep{Almgren2010}.  

That neutrino transport is so integral to the core-collapse supernova explosion mechanism, and yet so overwhelming in its computational demands, has motivated our choice: GenASiS\ is being developed with an eye towards cell-by-cell AMR \citep[e.g.][]{Khokhlov1998}. 
Block-structured approaches \citep{Berger1984,Berger1989,MacNeice2000} are more common, but are not necessarily the best choice if neutrino transport is a major focus. 
One initial motivation for block-structured AMR is the ability to deploy existing `unigrid' time-explicit hydrodynamics solvers on individual regular cell blocks. 
Moreover, there are some efficiencies associated with the use of predictable basic building blocks.
In a multiphysics context, however, cell-by-cell AMR may have some advantages.  
It is not conceptualized around local time-explicit solvers, so that it is arguably more amenable to elliptic and other global solvers, including those that are time-implicit.
It is of course possible to develop global solvers in block-structured AMR---for instance for gravity \citep[e.g.][]{OShea2005,Ricker2008,Almgren2010} and radiation \citep[e.g.][]{Rijkhorst2006,Wise2011,Zhang2011}, but it is not clear that this is the most natural environment imaginable for them.
Moreover, the computational cost per cell becomes very high if one aims (at least eventually, if not initially) towards the high dimensionality of the full neutrino phase space (position space plus momentum space), and towards very large nuclear reaction networks. 
In this case cell-by-cell AMR has a potentially significant advantage in requiring a smaller number of total cells for a given accuracy, in part because of the more fine-grained control over cell division and placement, and also because the use of many blocks at a given level of refinement in block-structured AMR leads to a larger ratio of ghost to computational cells.  

While a handful of other codes used in astrophysics or cosmology do use cell-by-cell AMR \citep[e.g.][]{Khokhlov1998,Teyssier2002,Gittings2008}, our approach in GenASiS\ has at least one notable difference. 
In arranging storage and writing solvers, rather than addressing the oct-tree as a single mesh consisting of the union of leaf cells at all levels of refinement, GenASiS\ has a more explicit level-by-level orientation. 
This facilitates multigrid approaches to elliptic and time-implicit global solves. 
Writing solves for single levels, with interactions between levels handled separately, restores some of the flavor of simplicity of the independent solves on individual blocks featured in block-structured AMR; at the same time, treating entire levels at once eliminates the drawbacks of having to stitch together results from multiple blocks at the same level.  

The high-level structure of the core of GenASiS\ is sketched in Figure~\ref{fig:GenASiS_Structure}. Solid lines outline relationships in the source code directory hierarchy, and dashed arrows indicate compilation dependencies. 
{\tt Modules} and {\tt Programs} are the two highest-level divisions shown.
{\tt Modules} comprises the classes forming the central functionality of GenASiS. Each class is defined in a single Fortran {\tt module}. 
Fortran {\tt program}s---drivers that invoke these classes, and begin the execution of some particular computational task---are collected under {\tt Programs}.  
Both {\tt Modules} and {\tt Programs} contain divisions labeled {\tt Physics}, {\tt Mathematics}, and {\tt Basics}. 
On the {\tt Modules} side, these categories contain class implementations; dual to these on {\tt Programs} side, in almost complete one-to-one correspondence, are {\em unit tests} that exercise the capabilities and provide example usage of the individual classes.
As indicated by the dashed dependency arrows, the unit tests first depend upon their corresponding classes: {\tt Basics} unit tests depend on {\tt Basics} classes; {\tt Mathematics} unit tests depend on {\tt Mathematics} classes, and through these also on {\tt Basics} classes; and so on. 
{\tt Programs} also has another division---{\tt Applications}---intended for drivers aimed at purposes beyond unit tests, ranging from the solution of simple physical test problems to the execution of production-scale multiphysics research simulations.
These ultimately depend on all {\tt Modules}---{\tt Physics}, {\tt Mathematics}, and {\tt Basics}.

With the top-down overview afforded by Figure~\ref{fig:GenASiS_Structure} in mind, this paper begins our series on GenASiS\ by describing the fundamental capabilities that underpin implementations of the solvers and physics needed for the simulation of core-collapse supernovae and other astrophysical systems. 
In particular we discuss the major code division {\tt Basics} and one of the major divisions of {\tt Mathematics}: cell-by-cell refinable {\tt Manifolds}.
Two sample problems illustrate our object-oriented approach and exercise the capabilities of the {\tt Basics} and {\tt Manifolds} divisions of GenASiS, including an example with a fixed multilevel grid of a type suitable for gravitational collapse. 

\section{BASICS}
\label{sec:Basics}

The {\tt Basics} division of GenASiS\ contains some utilitarian functionality for large-scale simulations on distributed-memory supercomputers. 
(For the place of {\tt Basics} in the overall scheme of GenASiS, see Figure \ref{fig:GenASiS_Structure}.) 
Its content, as illustrated in the left diagram of Figure~\ref{fig:Basics_Structure}, includes the divisions {\tt VariableManagement}, {\tt Display}, {\tt MessagePassing}, {\tt FileSystem}, and {\tt Runtime}. 
Two of these---{\tt Display} and {\tt Runtime}---are framed with boxes of thinner linewidth;
these are `leaf' divisions in the sense that they contain no further subdirectories, but only individual files defining classes.
The right diagram in Figure \ref{fig:Basics_Structure} shows the structure within {\tt VariableManagement}, which includes the leaf divisions {\tt Specifiers}, {\tt ArrayOperations}, {\tt ArrayArrays}, and {\tt VariableGroups}.
The middle right diagram in Figure \ref{fig:Basics_Structure} shows the structure within {\tt MessagePassing}, which includes the leaf divisions {\tt MessagePassingBasics}, {\tt PointToPoint}, and {\tt Collective}.
The middle left diagram in Figure \ref{fig:Basics_Structure} shows the structure within {\tt FileSystem}, which includes the leaf divisions {\tt FileSystemBasics}, {\tt GridImageBasics}, {\tt CurveImages}, {\tt StructuredGridImages}, and {\tt UnstructuredGridImages}.
In these illustrations the compilation order is from bottom to top; thus the dependencies essentially flow in reverse, from top to bottom.

Before proceeding further we say a bit more about the way in which GenASiS's classes are arranged and accessed.
In discussing the organization of GenASiS, sketched in Figures \ref{fig:GenASiS_Structure}, \ref{fig:Basics_Structure}, and similar figures to follow, we refer to the boxed entities as `divisions' of the code. 
They are not classes, but subdirectories in the source code tree; thus they comprise and compose hierarchically arranged groups of classes. 
However, even though these divisions are not classes per se, we nevertheless mirror their hierarchical structure in Fortran {\tt module}s that make groups of classes more convenient to access, as follows.
In each directory there is a file defining a Fortran {\tt module} whose name corresponds to the division or directory, and which contains all the classes (for a leaf division) or groups of classes (for a non-leaf division) within that directory. 
For instance, in the {\tt Specifiers} directory (see the right diagram in Figure \ref{fig:Basics_Structure}) there is a file {\tt Specifiers.f90} defining a {\tt module Specifiers} containing Fortran {\tt use} statements invoking each individual class in that directory.
Then the classes in the {\tt ArrayOperations} directory (the next to be compiled) need only the statement {\tt use Specifiers} to access all the classes in that grouping.
Similarly, the directory {\tt VariableManagement} contains a file {\tt VariableManagement.f90} as follows:
\begin{lstlisting}[frame=tb,numbers=left,numbersep=5pt,xleftmargin=10pt]
module VariableManagement
  use Specifiers
  use ArrayOperations
  use ArrayArrays
  use VariableGroups
end module VariableManagement
\end{lstlisting}
Then, in the classes grouped under, for example, {\tt Display} (see the left diagram in Figure \ref{fig:Basics_Structure}), the statement {\tt use VariableManagement} suffices to access everything contained within that division.

In this section we illustrate some of the functionality available through {\tt Basics}---and also give a flavor of our use of the object-oriented approach---in the context of a simple example problem.


\subsection{A Simple Problem and an Object-Oriented Solution}
\label{subsec:SimpleProblem}

The following tasks define a sample problem we use to give a flavor of our object-oriented approach and illustrate the functionality available through {\tt Basics}: 
(a) divide a one-dimensional periodic lattice of equidistant points into segments assigned to the several processes of a distributed-memory parallel program; 
(b) initialize a sine wave whose periodicity is consistent with that of the lattice; 
(c) exchange `ghost' data points between adjacent lattice segments assigned to different processes;
(d) compute the derivative of the sine wave with the familiar second-order centered finite difference formula; 
(e) evaluate the sine wave and its derivative at a test point;
(f) as if for a parameterized phase space type of plot, form coordinate pairs of the sine wave and its derivative for each point;
and (g) write output to the file system in a format suitable for visualization.
We have contrived a Fortran {\tt program} (Listing~\ref{lst:Program_SWD_Outline}) and two supporting {\tt module}s (Listings~\ref{lst:Module_SWD_Outline} and \ref{lst:Module_WD_Outline}) to accomplish these tasks using the GenASiS\ {\tt Basics} classes. 
(Obviously this is not a serious application, and for such a simple problem, what follows will seem---like killing a fly with a sledgehammer---as rather more than necessary. 
But overkill on a simple example now will make later discussions of more sophisticated parts of GenASiS\ more tractable by illustrating some of the object-oriented mechanisms deployed to more proportionate purpose as the functionality and applicability we discuss widen.)
We do not show every detail, but only selected snippets that give a few examples and illustrate some points we wish to make about our object-oriented approach and, in Section \ref{subsec:BasicsFunctionality}, the functionality available in GenASiS\ {\tt Basics} classes.

The {\tt program} outlined in Listing~\ref{lst:Program_SWD_Outline} accomplishes the tasks required by our simple problem.
\begin{lstlisting}[float,frame=tb,numbers=left,numbersep=5pt,xleftmargin=10pt,label=lst:Program_SWD_Outline,caption={Outline of {\tt program SineWaveDerivative}}]
program SineWaveDerivative

  use Basics
  use SineWaveDerivative_Form  !-- See Listing \ref{lst:Module_SWD_Outline}
  implicit none

  type ( SineWaveDerivativeForm ) :: SWD  !-- See Listing \ref{lst:Module_SWD_Outline}
  
  !-- Program initialization omitted; see Section~\ref{subsubsec:Runtime}

  call SWD
  call SWD
  call SWD
  call SWD
  call SWD
  call SWD

  !-- Program finalization omitted

end program SineWaveDerivative
\end{lstlisting}
The statement {\tt  use Basics} in line~3 gives access to all the classes in the code divisions pictured in Figure~\ref{fig:Basics_Structure}.
In line~7 we declare an object {\tt SWD}\footnote{The identifier {\tt SWD} is an acronym for {\tt SineWaveDerivative}. Such acronymic naming, especially for objects, is a very common practice of ours.}
of class {\tt SineWaveDerivativeForm},\footnote{Our convention is to use the suffix {\tt Form} in class names. 
We normally define one class per Fortran {\tt module}, with the suffix preceded by an underscore ({\tt \_Form}) in the {\tt module} name to avoid a collision with the corresponding derived {\tt type} name.}
whose definition is available through the {\tt use} statement in line~4.
The {\tt subroutine} calls in lines~11-16 are invocations of methods belonging to the object {\tt SWD}.\footnote{The syntax [variable] \% [procedure] to reference methods is new to Fortran 2003. 
Procedures---either {\tt subroutine}s or {\tt function}s---are attached to derived {\tt type}s through one of two new mechanisms: {\tt procedure} pointer components, or, as we shall see for instance in Listing \ref{lst:Module_SWD_Outline}, type-bound {\tt procedure} declarations \citep[e.g.][]{Reid2007}. 
By default the object that owns the method is the first---and so-called {\em passed}---argument to the invoked procedure, giving that procedure access to the object's members.
Additional arguments follow in parentheses. In lines~11-16 the parentheses following the method names are empty because the object---in this case {\tt SWD}---is the only argument to those {\tt subroutine}s.}  
The comments in these lines correlate the tasks outlined in the previous paragraph with the method invocations.
  
Evidently we must consider the class {\tt SineWaveDerivativeForm} defined in the {\tt module} sketched in Listing~\ref{lst:Module_SWD_Outline} to begin to see how these tasks are implemented.
\begin{lstlisting}[float,frame=tb,numbers=left,numbersep=5pt,xleftmargin=10pt,label=lst:Module_SWD_Outline,caption={Outline of {\tt module SineWaveDerivative\_Form}. Used at line~4 of Listing~\ref{lst:Program_SWD_Outline}.}]
module SineWaveDerivative_Form

  use Basics
  use WaveDerivative_Template  !-- See Listing \ref{lst:Module_WD_Outline}
  implicit none
  private

  type, public, extends ( WaveDerivativeTemplate ) :: SineWaveDerivativeForm  
    type ( MeasuredValueForm ) :: Amplitude
  contains
    procedure, private, pass :: Initialize_SWD  !-- See Listing \ref{lst:Subroutine_Initialize_SWD}
    generic, public :: Initialize => Initialize_SWD
    procedure, private, pass :: Waveform
  end type SineWaveDerivativeForm

contains

  !-- Definitions of {\tt subroutine Initialize\_SWD} and {\tt function Waveform} omitted

end module SineWaveDerivative_Form
\end{lstlisting}
One member of {\tt SineWaveDerivativeForm} is declared in line 9, the {\tt Amplitude} of the sine wave. 
It is an object of class {\tt MeasuredValueForm}; as we will discuss in Section \ref{subsec:BasicsFunctionality}, this allows {\tt Amplitude} to have a unit of measure associated with it.
Lines 11-13 are declarations related to two of the methods of {\tt SineWaveDerivativeForm}.\footnote{The procedures named in lines 11 and 13 have the {\tt  private} attribute, which means they are not accessible outside this {\tt module}. The default, which we usually affirm explicitly, is {\tt public} availability. 
These methods also have the {\tt pass} attribute, which affirms the default behavior that the object owning the method is passed as the first argument \citep[e.g.][]{Reid2007}.}
The reason for the {\tt generic} statement in line 12, which provides a {\tt public}ly available alias {\tt Initialize} for the {\tt private} method {\tt Initialize\_SWD} declared in line 11, will be made clear shortly.

Noting that of the methods invoked in lines 11-16 of Listing~\ref{lst:Program_SWD_Outline}, only one of them---{\tt Initialize}---is made publicly available (via line~12) in Listing~\ref{lst:Module_SWD_Outline}, we turn our attention to the indication in line~8 of Listing~\ref{lst:Module_SWD_Outline} that this class {\tt extends} the class {\tt WaveDerivativeTemplate} outlined in Listing~\ref{lst:Module_WD_Outline}. 
\begin{lstlisting}[float,frame=tb,numbers=left,numbersep=5pt,xleftmargin=10pt,label=lst:Module_WD_Outline,caption={Outline of {\tt module WaveDerivative\_Template}. Used at line~4 of Listing~\ref{lst:Module_SWD_Outline}.}]
module WaveDerivative_Template

  use Basics
  implicit none
  private

  type, public, abstract :: WaveDerivativeTemplate
    !-- Declaration of members omitted
  contains
    procedure, private, pass :: Initialize_WD     !-- See Listing \ref{lst:Subroutine_Initialize_WD}
    generic, public :: Initialize => Initialize_WD
    procedure, public, pass :: ExchangeVariables  !-- See Listing \ref{lst:Subroutine_ExchangeVariables}
    procedure, public, pass :: ComputeDerivative
    procedure, public, pass :: EvaluateVariables  !-- See Sections~\ref{subsubsec:VariableManagement}, \ref{subsubsec:Display}, \ref{subsubsec:MessagePassing}
    procedure, public, pass :: PairVariables      !-- See Listing \ref{lst:Subroutine_PairVariables}, and Section~\ref{subsubsec:VariableManagement}
    procedure, public, pass :: WriteVariables     !-- See Section~\ref{subsubsec:VariableManagement}
    procedure ( WaveformInterface ), private, pass, deferred :: Waveform
  end type WaveDerivativeTemplate

  !-- Declaration of the {\tt interface} for {\tt function WaveformInterface} omitted
  
contains

  !-- Definitions of {\tt subroutine}s {\tt Initialize\_WD}, {\tt ExchangeVariables},  {\tt ComputeDerivative}, 
  ! \  \   {\tt EvaluateVariables}, {\tt PairVariables}, and {\tt WriteVariables} omitted

end module WaveDerivative_Template
\end{lstlisting}
The extension of {\tt WaveDerivativeTemplate} to {\tt SineWaveDerivativeForm} is an example of the object-oriented principle of inheritance. 
The idea in this particular example is that {\tt WaveDerivativeTemplate} contains members and methods useful to solve this simple example problem for any waveform, while {\tt SineWaveDerivativeForm} adds only the elements necessary to specify that the waveform be a sine wave.
{\tt SineWaveDerivativeForm} is said to inherit all the members suggested by line~8 of Listing~\ref{lst:Module_WD_Outline}, and all the methods declared in lines~10-17.
This inheritance is what allows, for instance, the object {\tt SWD} of class {\tt SineWaveDerivativeForm} declared in line~8 of the {\tt program} in Listing {\ref{lst:Program_SWD_Outline} to make the method invocations given in lines~11-16 of that listing.

A peculiarity of this particular example is that {\tt WaveDerivativeTemplate} is {\tt abstract} (as specified in line~7 of Listing~\ref{lst:Module_WD_Outline}), as required by the fact that its method {\tt Waveform} is {\tt deferred} (as declared in line~17 of that listing).\footnote{Our convention is to use the suffix {\tt Template} for an {\tt abstract} class instead of the suffix {\tt Form} we use for typical classes.}
This means that, unlike all the other methods declared in lines~10-17 of Listing~\ref{lst:Module_WD_Outline}---which are defined or `fleshed out' in this class definition, as suggested by lines~24-25---the method {\tt Waveform} is only given an {\tt interface}, or specification of its argument list, as suggested by line~20.
In this example the method {\tt Waveform} is not fleshed out until {\tt WaveDerivativeTemplate} is extended to {\tt SineWaveDerivativeForm} in Listing~\ref{lst:Module_SWD_Outline}, in which the previous {\tt deferred} declaration is said to be {\em overridden} by the new declaration in line~13, and the method is defined as suggested in line~18.
We emphasize that inheritance by type extension does not only occur in going from an {\tt abstract} class to a non-{\tt abstract} one, but between any classes at all---and indeed through any number of sequential extensions.

Finally, we note multiple ways in which polymorphism (use of the same name to refer to different variations of the same basic concept) are exhibited in this example.\footnote{Amusingly, the object-oriented principle of polymorphism is itself polymorphic, manifesting itself in multiple ways!}
Consider first the {\em overloading} of the method {\tt Initialize}, which by virtue of the {\tt generic} statements in line 12 of Listing~\ref{lst:Module_SWD_Outline} and line 11 of Listing~\ref{lst:Module_WD_Outline}, ends up as an alias for two different {\tt subroutine}s named {\tt Initialize\_WD} and {\tt Initialize\_SWD}. 
This is possible because these two routines have different argument lists, which allows the system to {\em resolve} any particular call (i.e. route it to the correct {\tt subroutine}) by finding the aliased routine with the appropriate argument list.
To see how this works, and to then encounter additional manifestations of polymorphism, consider the {\tt subroutine}s {\tt Initialize\_SWD} and {\tt Initialize\_WD} in Listings~\ref{lst:Subroutine_Initialize_SWD} and \ref{lst:Subroutine_Initialize_WD}, which would fall in line~18 of Listing~\ref{lst:Module_SWD_Outline} of {\tt SineWaveDerivativeForm} and lines~24-25 of Listing~\ref{lst:Module_WD_Outline} of {\tt WaveDerivativeTemplate} respectively. 
\begin{lstlisting}[float,frame=tb,numbers=left,numbersep=5pt,xleftmargin=10pt,label=lst:Subroutine_Initialize_SWD,caption={{\tt subroutine Initialize\_SWD.} Fits in line~18 of Listing~\ref{lst:Module_SWD_Outline}. Called at line~12 of Listing~\ref{lst:Program_SWD_Outline}.}]
subroutine Initialize_SWD ( SWD )

  class ( SineWaveDerivativeForm ), intent ( inout ) :: SWD
  type ( UnitForm ) :: SineWaveUnit

  SineWaveUnit    = UNIT
  SWD
  !-- Statement overriding this default with a value read from parameter file or command line omitted

  call SWD

end subroutine Initialize_SWD
\end{lstlisting}
\begin{lstlisting}[float,frame=tb,numbers=left,numbersep=5pt,xleftmargin=10pt,label=lst:Subroutine_Initialize_WD,caption={{\tt subroutine Initialize\_WD}. Fits in lines~24-25 of Listing~\ref{lst:Module_WD_Outline}. Called at line~10 of Listing~\ref{lst:Subroutine_Initialize_SWD}.}]
subroutine Initialize_WD ( WD, WaveUnit )

  class ( WaveDerivativeTemplate ), intent ( inout ) :: WD
  type ( UnitForm ), intent ( in ) :: WaveUnit

  !-- Declaration of other variables, and other statements---some of which are in Listings \ref{lst:CoordinateInitialization} and \ref{lst:VariableGroupInitialization}---omitted
  
  !-- {\tt X} and {\tt W} are rank-one {\tt real} arrays storing the lattice point coordinates and wave values on those points respectively 
  W = WD
  
end subroutine Initialize_WD
\end{lstlisting}
Notice that {\tt Initialize\_SWD} has only one argument while {\tt Initialize\_WD} has two. 
This difference allows the call in line~11 of the {\tt program} in Listing~\ref{lst:Program_SWD_Outline} to be resolved to {\tt Initialize\_SWD} (with the object {\tt SWD} to the left of the {\tt \%} symbol passed as the only argument), and the call in line~10 of Listing~\ref{lst:Subroutine_Initialize_SWD} to be resolved to {\tt Initialize\_WD} (with {\tt SWD} as the first argument and {\tt SineWaveUnit} as the second argument).
The declaration of the argument {\tt WD} in line~3 of Listing \ref{lst:Subroutine_Initialize_WD} with the {\tt class} keyword instead of the {\tt type} keyword makes {\tt WD} a polymorphic variable, which means that this argument can be any extension (or extension of extension, and so on) of {\tt WaveDerivativeTemplate}.\footnote{That the {\tt class} keyword denotes a polymorphic variable in Fortran 2003 must not be confused with our use of the term `class' to refer to a construct that bundles together members and methods.} 
This is what allows a {\tt SineWaveDerivativeForm} object---ultimately, {\tt SWD} declared in line~7 of Listing {\ref{lst:Program_SWD_Outline}---to be the first argument of {\tt Initialize\_WD}.

Now for an important point regarding the utility of the object-oriented principles of inheritance and polymorphism for our purposes in GenASiS. 
The function call in line~9 of Listing~\ref{lst:Subroutine_Initialize_WD} of {\tt Initialize\_WD} (a method of {\tt WaveDerivativeTemplate} in Listing~\ref{lst:Module_WD_Outline}) sets the wave values on the coordinate points (see also the comment in line~8 of Listing~\ref{lst:Subroutine_Initialize_WD}).
As far as {\tt WaveDerivativeTemplate} is concerned this function could return anything; its definition is deferred, with only an interface specified (lines~17 and 20 of Listing \ref{lst:Module_WD_Outline}).
But when the argument {\tt WD} is an instance of {\tt SineWaveDerivativeForm}, the function call resolves to the overriding method specified by lines~13 and 18 of Listing~\ref{lst:Module_SWD_Outline} of {\tt SineWaveDerivativeForm}, returning values of a sine wave.
The point is that the principles of inheritance and polymorphism---exemplified here in the mechanisms of {\tt type} extension and method overriding---make it much easier to allow lower-level code (in this example, {\tt Initialize\_WD} in {\tt WaveDerivativeTemplate}) to access higher-level code (in this example, {\tt Waveform} in {\tt SineWaveDerivativeForm}).
This language functionality greatly facilitates, for instance, the separation of the {\tt Mathematics} and {\tt Physics} divisions of GenASiS. 
For example, solvers for generic classes of equations can be written in {\tt Mathematics}, and then invoked later by a range of different systems whose details are specified in {\tt Physics}.
This tremendously enhances the ease and transparency with which one can develop (for example) versatile and widely-applicable solvers.

(The class {\tt UnitForm} appearing in Listings~\ref{lst:Subroutine_Initialize_SWD} and \ref{lst:Subroutine_Initialize_WD}, and Lines~6 and 7 of Listing~\ref{lst:Subroutine_Initialize_SWD} which give a default setting of the sine wave amplitude, will be discussed shortly in Section \ref{subsec:BasicsFunctionality}.) 


\subsection{Illustrations of {\tt Basics} Functionality}
\label{subsec:BasicsFunctionality}

With this example program in mind we describe some of the functionality available in {\tt Basics}. In particular we discuss in bottom-up fashion (i.e. in order of compilation) the code divisions appearing in Figure \ref{fig:Basics_Structure}.

\subsubsection{VariableManagement}
\label{subsubsec:VariableManagement}

{\tt VariableManagement} is the first code division appearing in the diagram on the left side of Figure \ref{fig:Basics_Structure}. 
It is one of the {\tt Basics} divisions that is not a leaf division. 
Its subdivisions {\tt Specifiers}, {\tt ArrayOperations}, {\tt ArrayArrays}, and {\tt VariableGroups} are shown in the rightmost diagram of Figure \ref{fig:Basics_Structure}. 
We discuss each in turn.

\paragraph{Specifiers}

The first division of {\tt VariableManagement} is {\tt Specifiers}, which contains classes used in the specification of number and character variables. 
This includes some default parameters to be used as {\tt kind} and {\tt len} specifiers for Fortran intrinsic data types. 
It also includes mathematical and physical constants and a means of dealing with units. 

Consider for example the precision and range of numerical variables, which can be specified with a Fortran {\tt kind} parameter. 
Note for instance the declaration of one of the members of {\tt WaveDerivativeTemplate},\footnote{As a member of {\tt WaveDerivativeTemplate}, at for instance line~6 of Listing~\ref{lst:Subroutine_Initialize_WD} this variable would be referenced as {\tt WD \% Wavelength}.} which falls in line~8 in  Listing~\ref{lst:Module_WD_Outline}:
\begin{lstlisting}
  real ( KDR ) :: Wavelength
\end{lstlisting}
The parameter {\tt KDR} is a nickname, provided by {\tt Specifiers}, for {\tt KIND\_DEFAULT \% REAL}, one of the members of the object {\tt KIND\_DEFAULT} declared and initialized in {\tt Specifiers}.\footnote{By convention, when we write a class name that is not an acronym in all caps (for instance {\tt KIND\_DEFAULT}) it refers to a {\em singleton}, a special type of object intended to be the only instance of that class ever declared. 
When we write a class member (e.g. the {\tt REAL} in {\tt KIND\_DEFAULT \% REAL}) in all caps it denotes a parameter that is intended to remain unchanged throughout program execution.
All-caps class members remaining unchanged is at the mercy of convention, except in the case of a singleton object declared with the {\tt parameter} attribute, in which case all of its members automatically remain unchanged.}
Fortran {\tt kind} parameters also specify the precision and range of explicit numbers, as in the {\tt 1.0\_KDR} in line~7 of Listing~\ref{lst:Subroutine_Initialize_SWD}.
Normally {\tt KIND\_DEFAULT \% REAL} (and therefore its nickname {\tt KDR}) is set to correspond to double precision; but if one wanted to execute a program in, say, single precision (and assuming all real variables were declared to be of kind {\tt KDR}), changing this single parameter would accomplish the task.

{\tt Specifiers} also provides means for dealing with units, primarily for initialization and output.
These means include the classes {\tt MeasuredValueForm} and {\tt UnitForm} and a singleton object {\tt UNIT}.
The class {\tt MeasuredValueForm} has a {\tt real} member {\tt Number} and a {\tt character} member {\tt Unit} that specifies a unit of measure (e.g. {\tt 'm'} for meters).\footnote{Our apologies for using `unit' in three ways here: first, as the name of a {\tt character} member of {\tt MeasuredValueForm}; second, in the name of the class {\tt UnitForm}; and third, as the name of the singleton object {\tt UNIT}.
Which sense is meant in a particular expression should be clear in context.} 
The class {\tt UnitForm} is meant only for units of measure, and not for typical measured quantities. 
It is an extension of {\tt MeasuredValueForm} that has an additional {\tt character} member {\tt Label}, the need for which will be explained momentarily.
The members of the singleton {\tt UNIT} are objects of {\tt UnitForm} for many units of measure. 
Now consider one of these members, {\tt UNIT \% CENTIMETER}, as an example of the need for the member {\tt Label} of {\tt UnitForm}.
The need for this extra member (relative to the members {\tt Number} and {\tt Unit} of the parent class {\tt MeasuredValueForm}) arises from that fact that in its internal computations on plain {\tt real} variables, GenASiS\ assumes all numbers are given in terms of a single fundamental unit. 
To be specific, a geometrized system of units is currently employed in which the meter is fundamental.
In terms of this fundamental unit, 1 cm = $10^{-2}$ m. 
The right-hand side of this expression is reflected in the {\tt MeasuredValueForm} members {\tt Number} and {\tt Unit} of {\tt UnitForm}, i.e. {\tt UNIT \% CENTIMETER \% Number = 1.0e-2\_KDR} and {\tt UNIT \% CENTIMETER \% Unit = 'm'}; 
while the {\tt Label} member is associated with the unit in question, i.e. {\tt UNIT \% CENTIMETER \% Label = 'cm'}.
The members {\tt Number} and {\tt Unit} of {\tt UnitForm} are used to convert a number into or out of GenASiS's fundamental unit on input or output, while the {\tt Label} member is used to label the output of a number in the given unit.\footnote{We emphasize that type extension via the Fortran 2003 {\tt extends} keyword (as in line~8 of Listing~\ref{lst:Module_SWD_Outline}) gives the extended type direct access to the members of the parent type. 
For example, if {\tt CoordinateUnit} is declared to be of type {\tt UnitForm}, we can refer directly to the members {\tt CoordinateUnit \% Number} and {\tt CoordinateUnit \% Unit} inherited from the parent type {\tt MeasuredValueForm}, in addition to the member {\tt CoordinateUnit \% Label} added to the extended type.} 

We are now in a position to see how initializations with units can be coded directly. (Examples of initialization with units when reading numbers from parameter files or as command line options, and of output with units to the screen or to a file, will be discussed in connection with the {\tt FileSystem}, {\tt Runtime}, and {\tt Display} portions of {\tt Basics}.)
Here is the declaration of two other members of {\tt WaveDerivativeTemplate}, which falls in line~8 in  Listing~\ref{lst:Module_WD_Outline}:
\begin{lstlisting}
    type ( MeasuredValueForm ) :: MinCoordinate, MaxCoordinate
\end{lstlisting}
Listing~\ref{lst:CoordinateInitialization}, which would fall at line~6 of Listing~\ref{lst:Subroutine_Initialize_WD} of {\tt subroutine Initialize\_WD}, sets the default endpoints and units of the periodic lattice, and uses these endpoints (together with an {\tt integer} member {\tt nWavelengths} specifying the number of wavelengths of the periodic waveform to set on the lattice) to set the wavelength.
\begin{lstlisting}[float,frame=tb,numbers=left,numbersep=5pt,xleftmargin=10pt,label=lst:CoordinateInitialization,caption={Coordinate and wavelength initialization in {\tt Initialize\_WD}. Fits in line~6 of Listing~\ref{lst:Subroutine_Initialize_WD}.}]
  type ( UnitForm ) :: CoordinateUnit
  !-- Other declarations and statements omitted

  CoordinateUnit     = UNIT
  WD
  WD
  !-- Statement overriding these defaults with values read from a parameter file or command line omitted
  !-- Display to screen omitted; see Section~\ref{subsubsec:Display}

  WD
\end{lstlisting}
The arithmetic and assignment operations involving derived {\tt type}s in lines~5 and 6 of Listing~\ref{lst:CoordinateInitialization} behave in expected ways: the explicit {\tt real} number multiplies the {\tt Number} member on the right-hand side, and the {\tt Label} member in the {\tt UnitForm} on the right-hand side is simply ignored in the assignment to the {\tt MeasuredValueForm} on the left-hand side.\footnote{These and like results are achieved by overloading of the arithmetic and assignment operators.} 
Multiplication of an explicit number by one of the {\tt UNIT} members puts the coordinates of the lattice endpoints into GenASiS's fundamental working unit, so that one can assign their difference to the plain {\tt real} member {\tt Wavelength} in line~10. 
(In the assignment of a {\tt MeasuredValueForm} to a {\tt real}, the {\tt Unit} member is simply ignored.)
Another example is lines~6 and 7 of Listing~\ref{lst:Subroutine_Initialize_SWD}, in which the member {\tt Amplitude} of {\tt SineWaveDerivativeForm} (declared at line~9 of Listing~\ref{lst:Module_SWD_Outline}) receives a default initialization to a dimensionless value of unity.

Finally, we note that mathematical and physical constants can be accessed through the singleton {\tt CONSTANT}.
For instance, the number $\pi$ is referenced in the expression
\begin{lstlisting}
2.0_KDR * CONSTANT
\end{lstlisting}
for the wavenumber in {\tt function Waveform} (whose definition is noted at line 18 of Listing~\ref{lst:Module_SWD_Outline} for {\tt SineWaveDerivative\_Form}). 


\paragraph{ArrayOperations}

The second division of {\tt VariableManagement} is {\tt ArrayOperations}, which includes some basic operations on arrays.
Some of these---such as {\tt Clear} and {\tt Copy}---are merely overloaded names for {\tt subroutine}s whose arguments are intrinsic data types stripped of any potentially complicating contexts or attributes that might inhibit compiler optimization.\footnote{The modules defining the these two operations are {\tt Clear\_Command} and {\tt Copy\_Command}. 
We use the suffix {\tt Command} instead of {\tt Form} in naming a module when it is more convenient to simply implement an action, as opposed to a class with both members and methods. 
Slavish adherence to object-oriented principles would implement such a case as a class with a (perhaps overloaded) method and no members.
In the case of pure actions this is more cumbersome than necessary, leading us to revert to a more procedural rather than object-oriented perspective for {\tt Command}s.
Therefore the mechanism of overloading in this case is not the Fortran 2003 {\tt generic} statement used in connection with type-bound procedures; instead it is the use, possible since Fortran 90, of an {\tt interface} block to define a generic name for a set of procedures.
Nevertheless the basic idea of the overloading is still the same, i.e. use of a single name to refer to several possible routines, distinguished by and resolved according to their distinct argument lists.} 
The purpose of these wrappers is to expose elementary variables to the compiler so that it might use fast methods (such as a system call to set values in memory) it might not otherwise employ if, say, the array has the {\tt pointer} attribute, is deeply embedded in a complicated expression involving a sequence of derived {\tt type} member references, etc.
Other array operations have a more substantial purpose, such as {\tt Sort} and {\tt Search}.\footnote{These are defined in the modules {\tt Sort\_Command} and {\tt Search\_Command}.} 
The {\tt Search} operation, for instance, is used in task (e) of our sample problem (evaluate the sine wave and its derivative at a test point).
The method {\tt EvaluateVariables} of {\tt WaveDerivativeTemplate} includes the line
\begin{lstlisting}
  call Search ( Coordinate, X_Evaluate, iE )
\end{lstlisting}
The first two arguments are inputs: {\tt Coordinate} is a rank-one {\tt real} array containing the coordinates of the lattice points, and {\tt X\_Evaluate} is the {\tt real} coordinate of the test point at which the sine wave and its derivative are to be evaluated. 
Upon return, the last argument {\tt iE} is set to the index of {\tt Coordinate} for which {\tt Coordinate ( iE ) <= }{\tt X\_Evaluate < }{\tt Coordinate ( iE + 1 )}. 
The search having bracketed the test point between array indices {\tt iE} and {\tt iE + 1}, the coordinate and sine wave and derivative values at these points can then be used to perform the desired evaluations by linear interpolation.

\paragraph{ArrayArrays}

The next division of {\tt VariableManagement} is {\tt ArrayArrays}, which contains classes that can be used to form arrays of arrays.
One application is the construction of so-called {\em ragged arrays}. 
In the rank-two case a ragged array can be conceptualized (following the Fortran row-major perspective) as a matrix whose columns have different numbers of elements.
This can be useful when one wants an indexed collection of rank-one arrays that are not all the same length.
Even if they are all the same length, one might want to conceptualize the data primarily as a collection of rank-one objects rather than as a rank-two array.
This is the way our sample program handles task (f)---set up coordinate pairs of the sine wave and its derivative at each point, as if for a phase space type of plot.
The method {\tt PairVariables} of {\tt WaveDerivativeTemplate} is outlined in Listing~\ref{lst:Subroutine_PairVariables}.
\begin{lstlisting}[float,frame=tb,numbers=left,numbersep=5pt,xleftmargin=10pt,label=lst:Subroutine_PairVariables,caption={{\tt subroutine PairVariables}. Fits in lines~24-25 of Listing~\ref{lst:Module_WD_Outline}.}]
subroutine PairVariables ( WD )

  class ( WaveDerivativeTemplate ), intent ( inout ) :: WD

  integer ( KDI ) :: iP  !-- iPair
  !-- Declaration related to output to screen omitted
  type ( ArrayReal_1D_Form ), dimension ( : ), allocatable :: VariablePair
  
  !-- Coordinate, Wave, and Derivative are rank-one arrays
  allocate ( VariablePair ( size ( Coordinate ) ) )
  do iP = 1, size ( VariablePair )
    call VariablePair ( iP )
    VariablePair ( iP )
    !-- Output to screen omitted
  end do

end subroutine PairVariables
\end{lstlisting}
We want one variable pair---a value of the wave paired with a value of the derivative---for each coordinate point.
Accordingly we declare a rank-one array {\tt VariablePair} of {\tt ArrayReal\_1D\_Form} (line 7) and allocate it to have the same number of elements as the array of coordinate values (line 10). 
The class {\tt ArrayReal\_1D\_Form} has an allocatable rank-one {\tt real} array member {\tt Value}. 
In line 12, calling the {\tt Initialize} method of each element of the array {\tt VariablePair} with the argument {\tt 2} causes the member {\tt Value} to be allocated to be two elements long.\footnote{We emphasize that the method {\tt Initialize} of {\tt ArrayReal\_1D\_Form} invoked in line 12 is completely distinct and separate from the {\tt Initialize} method of {\tt WaveDerivativeTemplate} and its extension {\tt SineWaveDerivativeForm} we have seen previously.}
Then the wave and its derivative at a single point are assigned to the member {\tt Value} of one of the {\tt VariablePair} elements in line 13.
Thus we have formed a rank-one array (one for each {\tt Coordinate} point; lines 7, 10) of rank-one arrays of length 2 (for wave/derivative coordinate pairs; lines 12-13).
As an aside, we note that the method {\tt Initialize} of {\tt ArrayReal\_1D\_Form} is overloaded; possible alternative arguments include an array shape (used here), an explicit array to be copied, or another {\tt ArrayReal\_1D\_Form} to be copied.
For instance, the single line
\begin{lstlisting}
  call VariablePair ( iP )
\end{lstlisting}
could have been used instead of lines 16-17.
In addition to {\tt ArrayReal\_1D\_Form} there are of course other analogous classes for arrays of other types and ranks.

\paragraph{VariableGroups}

The final division of {\tt VariableManagement} is {\tt VariableGroups}, which contains some infrastructure we use extensively in handling collections of variables, especially sets of related physical fields. 
In particular, the class {\tt VariableGroupForm} includes both metadata about the variables (names, units, etc.) and storage for the variable data itself. 
The data storage is a rank-two {\tt real} array member {\tt Value} of {\tt VariableGroupForm}; the first dimension indexes separate points, and the second dimension indexes individual variables. 

Our sample program uses an instance of  {\tt VariableGroupForm} to store the variables on its lattice points.
To prepare to understand its initialization we note the declarations of some other members of {\tt WaveDerivativeTemplate}, which fall in line~8 in  Listing~\ref{lst:Module_WD_Outline}:
\begin{lstlisting}
  integer ( KDI ) :: COORDINATE = 1, WAVE = 2, DERIVATIVE = 3, N_VARIABLES = 3, nPointsLocal
  type ( VariableGroupForm ) :: VariableGroup
\end{lstlisting}
Member {\tt VariableGroup} will store three variables, hence {\tt N\_VARIABLES = 3}. 
To avoid the confusing appearance of plain numerical indices, arrays over these variables will indexed by the members {\tt COORDINATE}, {\tt WAVE}, and {\tt DERIVATIVE}, as will be made more clear shortly.
The integer {\tt nPointsLocal} stores the number of local lattice points assigned to each segment of the distributed lattice (recall that the lattice as a whole is divided among the processes of a distributed-memory program). 
Listing~\ref{lst:VariableGroupInitialization}, which would fall at line~6 of Listing~\ref{lst:Subroutine_Initialize_WD} of {\tt subroutine Initialize\_WD}, shows the initialization of the {\tt VariableGroup} member.
\begin{lstlisting}[float,frame=tb,numbers=left,numbersep=5pt,xleftmargin=10pt,label=lst:VariableGroupInitialization,caption={Variable group initialization in {\tt Initialize\_WD}. Fits in line~6 of Listing~\ref{lst:Subroutine_Initialize_WD}.}]
  type ( UnitForm ), dimension ( WD
  character ( LDL ), dimension ( WD
  !-- Other declarations and statements omitted

  VariableUnit ( WD
  VariableUnit ( WD
  VariableUnit ( WD

  VariableName ( WD
  VariableName ( WD
  VariableName ( WD

  call WD
         ( [ WD
           VariableOption = VariableName, UnitOption = VariableUnit, ClearOption = .true. )
\end{lstlisting}
Lines~1 and 2 of Listing~\ref{lst:VariableGroupInitialization} initialize local variables used to prepare the variable metadata, i.e. the units and names of the three variables to be stored in member {\tt VariableGroup}. 
This metadata is then prepared in lines 5-7 and 9-11.
(We have already seen the declaration of {\tt WaveUnit} and {\tt CoordinateUnit} as instances of {\tt UnitForm} in previous listings.\footnote{Note the division of one instance of {\tt UnitForm} by another in the right-hand side of line 7.
We have overloaded the arithmetic operators in such a way that arithmetic can be performed on instances of {\tt MeasuredValueForm} and {\tt UnitForm}, including some basic handling of the consequences for their {\tt character} members expressing units of measure.})
The initialization of the {\tt VariableGroup} member then comes in the call in line 13 to the {\tt Initialize} method of {\tt VariableGroupForm}. 
The first argument is an array with two elements specifying the shape of the data storage: the first element gives the number of points (the number of points in the local lattice segment plus two ghost points), and the second element specifies the number of variables.
The optional arguments that follow specify the aforementioned variable metadata, and instruct that the data array should be `cleared' (initialized to 0.0). 

It is often convenient to use aliases to access the variable data in an instance of {\tt VariableGroupForm}. 
Consider for instance the rank-one arrays {\tt Coordinate}, {\tt Wave}, and {\tt Derivative} referred to in Listing~\ref{lst:Subroutine_PairVariables} of the method {\tt PairVariables} of {\tt WaveDerivativeTemplate}.
The following Fortran 2003 (nested) {\tt associate} constructs\footnote{This aliasing mechanism, new to Fortran 2003, avoids the need to declare pointer variables to be used as shorthand names and allows the compiler to produce more efficient code.} placed at line~9 in Listing~\ref{lst:Subroutine_PairVariables} arranges these aliases for use as rank-one arrays:
\begin{lstlisting}
  associate  ( VG => WD
  associate &
    ( Coordinate => VG
      Wave       => VG
      Derivative => VG
\end{lstlisting}
An appropriate number of {\tt end associate} statements should then follow line~15.

A {\tt VariableGroupForm} can have its own storage, as in the above example; or point to existing storage, perhaps for the purpose of selecting only a subset of variables from an existing variable group. 
For instance, in method {\tt WriteVariables} of {\tt WaveDerivativeTemplate} the two dependent variables (the wave and its derivative) are treated differently from the independent variable (the coordinate of the lattice points). 
Thus in {\tt WriteVariables} a new instance of {\tt VariableGroupForm} is set up as a local variable:
\begin{lstlisting}
  type ( VariableGroupForm ) :: VG_Curves
    
  call VG_Curves
    ( WD
\end{lstlisting}
Again, overloading of the {\tt Initialize} method allows variation in behavior. 
In this case the local object {\tt VG\_Curves} uses the member {\tt VariableGroup} of an instance of {\tt WaveDerivativeTemplate} as the basis of its initialization, and in fact will point to its existing data storage rather than allocate its own and perform a copy operation. 
But as the second argument suggests, only our dependent variables are `selected' in the new variable group {\tt VG\_Curves}.
Code written to use objects of {\tt VariableGroupForm}---such as the GenASiS\ I/O classes invoked by method {\tt WriteVariables} of our sample program---will then act only on the `selected' variables.

One may wonder why the trouble should be taken to use a class like {\tt VariableGroupForm} rather than work directly with explicitly declared arrays of variables.
The reason is that this simple construct greatly simplifies and renders more generic many tasks.
For instance, as just mentioned, the file system facilities described later are written in terms of variable groups, so that a user need not program such operations in detail for the many sets of variables that may arise in different problems.
Another example is the exchange of ghost cell data in multidimensional manifolds (Section \ref{sec:Manifolds}), which is also handled in terms of variable groups.
In writing versatile and reusable code it is much better to be able to simply set up a variable group rather than have to continually rewrite (for example) explicit ghost exchange or I/O library calls for the many different variables that may appear in, or be added to or removed from, various applications. 

Finally, we mention that the code division {\tt VariableGroups} also includes a class {\tt PackedVariableGroupForm} used to collect selected rows and columns of a {\tt VariableGroupForm} into contiguous memory for more efficient processing.

\subsubsection{Display}
\label{subsubsec:Display}

After {\tt VariableManagement} comes the leaf division {\tt Display} (see the diagram on the left side of Figure \ref{fig:Basics_Structure}), which contains infrastructure for displaying messages and variables to the standard output in a uniform and orderly way. 
The two central features available through {\tt Display} are the singleton object {\tt CONSOLE}, which has members that are used to control output to the standard display (`stdout'); 
and the heavily overloaded routine {\tt Show},\footnote{Defined in the module {\tt Show\_Command}.} which manages the hassles of displaying various data types with their respective formatting.

One issue that arises in a distributed-memory parallel computing environment is that, absent intervention by the programmer, output emerges from all processes in an asynchronous manner.
This typically results in unordered output of unwanted redundancy. 
The use of {\tt if} statements can avoid this by allowing only a particular process to output to the screen, but a proliferation of {\tt if} statements with every {\tt print} statement becomes tedious (and difficult to change for debugging purposes if a global variable specifying the process to be displayed is not used).
Thus {\tt CONSOLE} has a member {\tt DisplayRank} used by the {\tt Show} command to display output only from the process of the chosen rank, without the programmer having to continually supply {\tt if} statements.
If needed, the member {\tt DisplayRank} can be changed anywhere in the program.
Consider for instance task (e) in our sample program (evaluate the sine wave and its derivative at a test point), implemented by the method {\tt EvaluateVariables} of {\tt WaveDerivativeTemplate}. 
In general the test point will not fall in the lattice segment belonging to the default setting of {\tt DisplayRank}.
Assuming all processes are aware of {\tt EvaluateRank} (the rank of the process owning the relevant lattice segment), control of the display can be given to that process as follows:
\begin{lstlisting}
  DisplayRankOld = CONSOLE
  call CONSOLE
  !-- Perform function evaluations and display results
  call CONSOLE
\end{lstlisting} 
The second invocation of the {\tt SetDisplayRank} method returns control of the display to the process that previously had it.
Alternatively, instead of calling {\tt SetDisplayRank} one can supply an optional argument to the {\tt Show} command specifying, for that call only, which process should be visible.

Another issue related to output is that different pieces of information are of different levels of interest. 
In a production run one might want minimal output, while more verbose output may be useful for debugging.
{\tt CONSOLE} has a member {\tt Verbosity} that addresses this need.
It can be set to one of several levels denoted, in descending order of importance, by the {\tt CONSOLE} members {\tt ERROR}, {\tt WARNING}, {\tt INFO\_1}, {\tt ...,} {\tt INFO\_7}.
In a call to {\tt Show} one of these same members can be given as an optional argument; if the importance is thus tagged as less than {\tt Verbosity}, the output is suppressed. 
In this optional argument's absence an importance of {\tt INFO\_1} is assumed.

The {\tt Show} command provides a very simple interface for a wide range of formatted and labeled output.
There are provisions for intrinsic data types, objects of {\tt MeasuredValueForm}, and arrays of all these.
An object of {\tt UnitForm} can be supplied to modify the units (including the unit label itself) in which {\tt real} or {\tt MeasuredValueForm} numbers are displayed.
Consider for instance the following lines, which fit in line~8 of Listing~\ref{lst:CoordinateInitialization}:
\begin{lstlisting}
  call Show ( WD
  call Show ( WD
  call Show ( WD
  call Show ( WD
\end{lstlisting}
Recall that the {\tt MinCoordinate} and {\tt MaxCoordinate} members are of {\tt MeasuredValueForm}, and that their assignment in lines~5-6 of Listing~\ref{lst:CoordinateInitialization} puts them into GenASiS's internal geometrized units---in this case, meters.
This means that the first two {\tt Show} calls above will give output in this unit, tagged with the unit string {\tt 'm'}, and labeled with the strings {\tt 'MinCoordinate'} and {\tt 'MaxCoordinate'}. 
The output from the last two calls, which supply an additional argument {\tt CoordinateUnit} (set to {\tt UNIT \% CENTIMETER} in line~4 of of Listing~\ref{lst:CoordinateInitialization}), will be appropriately converted and tagged with the unit string {\tt 'cm'}.
The output will only come from the process specified by {\tt CONSOLE \% DisplayRank}, and the absence of the optional argument specifying importance means they will be regarded as {\tt INFO\_1}.
As another example consider the line
\begin{lstlisting}
  call Show ( Derivative, DerivativeUnit, 'Derivative', CONSOLE
\end{lstlisting}
Here {\tt Derivative} is an alias to the array of the wave derivative values. {\tt DerivativeUnit} is an alias to the value set in line~7 of Listing~\ref{lst:VariableGroupInitialization}; 
with the defaults of a dimensionless {\tt WaveUnit} (see Listing~\ref{lst:Subroutine_Initialize_SWD}) and {\tt CoordinateUnit} of centimeters in our sample program, {\tt DerivativeUnit} turns out to be inverse centimeters.
Thus this {\tt Show} call will output each value of the array of derivative values on a separate line, labeled with the unit string {\tt 'cm\textasciicircum-1'}.
The output will only appear from the process of rank {\tt CONSOLE \% DisplayRank}, and only if {\tt CONSOLE \% Verbosity} is high enough to include {\tt INFO\_3}.\footnote{The setting of {\tt CONSOLE \% Verbosity} for a program is done via a user-specified parameter file (see Section \ref{subsubsec:FileSystem}) and/or command-line arguments (see Section \ref{subsubsec:Runtime}).}


\subsubsection{MessagePassing}
\label{subsubsec:MessagePassing}

The next division of {\tt Basics} is {\tt MessagePassing} (see the diagram on the left side of Figure \ref{fig:Basics_Structure}), which contains some classes that abstract the data and methods useful for working in a message passing parallel computing environment.
As shown in the middle right diagram in Figure \ref{fig:Basics_Structure}, it has three subdivisions: {\tt MessagePassingBasics}, {\tt PointToPoint}, and {\tt Collective}.

Before describing the leaf divisions under {\tt MessagePassing} we note that one may ask why we bother to create and use message passing classes when MPI (Message Passing Interface) \citep{Gropp1999} is well-established, widely used, and seems likely to remain for many years the dominant library for message passing in distributed-memory applications. 
Our message passing classes do indeed use MPI `under the hood,' but in principle could be given other back ends without modifying their references in higher-level code.
Our interfaces are somewhat streamlined relative to direct MPI calls while providing methods for common use cases relevant to our system, and thus provide a simpler experience for the user. 
In terms of object-oriented design patterns \citep{Gamma1994}, our {\tt MessagePassing} classes provide a `fa\c{c}ade' that shields the user from concerns  such as the {\tt kind} and size of data types being sent and received (e.g. single precision vs. double precision, and beyond, in the case of {\tt real} data) and declaring many variables---handles, buffers, tags, and so on---for use in MPI calls. 
In our approach many of these needed variables are encapsulated in the classes we define. 
Finally, we note that the MPI library is mostly procedural in orientation and uses Fortran 77 semantics for its routines, sometimes requiring more arguments than necessary if the equivalent modern Fortran interfaces were available. 
Our interfaces relieve the user of these burdens.

\paragraph{MessagePassingBasics}
\label{subsubsec:MessagePassingBasics}

The most important class to mention here is {\tt CommunicatorForm}. 
Inspired (obviously) by the notion of an MPI communicator, this class connects a group of processes in a distributed-memory parallel program. 
The members of {\tt CommunicatorForm} include its {\tt Handle}, {\tt Size}, and {\tt Rank}. 
Its method {\tt Initialize} is overloaded to provide for either starting the message passing system, i.e. initializing MPI and thus creating the `global' or `world' communicator; or for creating a `subcommunicator' that connects only a subset of the processes connected by a parent communicator. 
The method {\tt Synchronize} pauses all processes until they reach the call to this method (i.e. an MPI `barrier'). 
For instance, recall the discussion in the second paragraph of Section \ref{subsubsec:Display} on changing the {\tt DisplayRank} during the execution of task (e) (evaluate the sine wave and its derivative at a test point). 
To ensure the output from two different processes remains ordered (and with the alias {\tt C} for the {\tt CommunicatorForm} object representing the global communicator), the line
\begin{lstlisting}
  call C
\end{lstlisting}  
is included before each of the two calls to the {\tt SetDisplayRank} method. 
Another method of {\tt CommunicatorForm} is {\tt Abort}, which provides an emergency exit from distributed-memory program execution.
Normal termination of the message passing system (i.e. finalization of MPI) occurs when the {\tt CommunicatorForm} object representing the global communicator goes out of existence.

\paragraph{PointToPoint} A single message from one process to another is a `point-to-point' communication. 
The division {\tt PointToPoint} provides classes for single incoming and outgoing messages of different intrinsic data types, as well as for arrays of such messages. 
These are used in task (c) of our sample program, as shown in Listing~\ref{lst:Subroutine_ExchangeVariables} for the method {\tt ExchangeVariables} of {\tt WaveDerivativeTemplate}.
\begin{lstlisting}[float,frame=tb,numbers=left,numbersep=5pt,xleftmargin=10pt,label=lst:Subroutine_ExchangeVariables,caption={{\tt subroutine ExchangeVariables}. Fits in lines~24-25 of Listing~\ref{lst:Module_WD_Outline}. Called at line~13 of Listing~\ref{lst:Program_SWD_Outline}.}]
subroutine ExchangeVariables ( WD )

  class ( WaveDerivativeTemplate ), intent ( inout ) :: WD

  integer ( KDI ) :: PreviousRank, NextRank
  type ( IncomingMessageArrayRealForm ) :: IMA
  type ( OutgoingMessageArrayRealForm ) :: OMA

  !-- Aliasing via associate construct omitted
  
  PreviousRank = modulo ( myRank - 1, nRanks )
  NextRank     = modulo ( myRank + 1, nRanks )

  call IMA
  call IMA

  call OMA
  OMA
  OMA
  call OMA

  call IMA
  Coordinate ( iIG ) = IMA
  Wave       ( iIG ) = IMA
  Coordinate ( iOG ) = IMA
  Wave       ( iOG ) = IMA

  call OMA
  
  !-- Display to screen omitted
  
end subroutine ExchangeVariables
\end{lstlisting}
The exchange of ghost point data requires each process to receive two messages and send two messages, one each to the processes owning the `previous' and `next' segments of the lattice (the ranks of these processes are set in lines 11 and 12).
Lines 6 and 7 declare objects {\tt IMA} and {\tt OMA} of classes {\tt IncomingMessageArrayRealForm} and {\tt OutgoingMessageArrayRealForm} respectively.
In this case the initializations of {\tt IMA} and {\tt OMA} (lines 14 and 17) are identical. 
The first argument is an object of {\tt CommunicatorForm}, here aliased as {\tt C}. 
That the subsequent arguments are arrays of length 2 indicates that {\tt IMA} and {\tt OMA} will each have two messages.
The second argument is an array of tags; in this case these are trivial because the fact that only one message will arrive from a given process means that there is no need to more specifically label (i.e. tag) them.
The third argument is an array of the ranks of the process ranks from (to) which the two messages are to be received (sent).
The fourth argument is an array containing the numbers of elements in each message; this is 2 for each message because we are sending two values for each ghost point, i.e. a coordinate value and a wave value.
Non-blocking `receives' of the incoming messages are posted by the single call in line 15.
Non-blocking `sends' are initiated (line 20) after the outgoing message buffers to the previous and next processes are loaded (lines 18 and 19) with data from the first (indexed by alias {\tt iF}) and last (indexed by alias {\tt iL}) points in the segment.
The {\tt WaitAll} method (line 22) ensures that the data have been received and are available in the incoming message buffers before they are stored in the inner (indexed by alias {\tt iIG}) and outer (indexed by alias {\tt iOG}) ghost points (lines 23-26).
A perfunctory wait for the outgoing messages to finish sending (line 28) completes the exchange.

\paragraph{Collective}

A communication that involves a group of processes is said to be `collective'. 
The division {\tt Collective} contains classes that provide for collective communications involving different intrinsic data types.
The collective communication class for a given data type is generic in the sense that its members are sufficient for broadcast, gather, scatter, and combined scatter-gather\footnote{`AllToAll' routines in MPI.} communications, and for collection operations (i.e. reductions of data); 
thus all of these are simply implemented as methods of a single class. 
For the gather and reduction methods, the absence or presence of the optional argument {\tt RootOption} in the initialization of a collective communication object determines whether or not the gathered or reduced value is distributed to all processes.\footnote{In MPI this corresponds to {\tt MPI\_Allgather} vs. {\tt MPI\_Gather}, or {\tt MPI\_Allreduce} vs. {\tt MPI\_Reduce}.}
This option should also be used in connection with broadcast and scatter calls: by definition, these communications involve a `root' process that is the source of data to be sent to all the processes in the group.
As an example of a collective communication in our sample program, the {\tt EvaluateVariables} method of {\tt WaveDerivativeTemplate} makes all processes aware, via a gather operation, of which process owns the test point and will perform the function evaluations (task (e); see also Section \ref{subsubsec:VariableManagement}, ``ArrayOperations''; and Section \ref{subsubsec:Display}). 
Consider the following lines:
\begin{lstlisting}
  call CO_Evaluate
  CO_Evaluate
  call CO_Evaluate
  EvaluateRank = maxloc ( CO_Evaluate
\end{lstlisting}
The object {\tt CO\_Evaluate} is an instance of the class {\tt CollectiveOperationIntegerForm}.
The first argument in its initialization in the first line is an object of {\tt CommunicatorForm}, aliased here as {\tt C}.
The second and third arguments in the initialization specify how large the outgoing and incoming buffers need to be.
In the second line, the integer {\tt myEvaluate} assigned to the outgoing buffer on each process has been set to {\tt 0} on all processes except the one owning the the test point, which sets it to {\tt 1}.
The gather operation is invoked in the third line; the absence of the optional argument {\tt RootOption} in the initialization means that the gathered value is distributed to all processes.
Then all processes can determine {\tt EvaluateRank} by finding the index of the incoming buffer element that is nonzero. 


\subsubsection{FileSystem}
\label{subsubsec:FileSystem}

As seen in the diagram on the left side of Figure \ref{fig:Basics_Structure}, the next division of {\tt Basics} is {\tt FileSystem}, whose classes handle I/O to disk.
The middle left diagram in Figure \ref{fig:Basics_Structure} shows divisions {\tt FileSystemBasics}, {\tt GridImageBasics}, {\tt CurveImages}, {\tt StructuredGridImages}, and {\tt UnstructuredGridImages}.  
The major classes in {\tt FileSystem} used to create and interact with files of various kinds have the term `stream' in their names, connoting a source or sink of data flow.
The names of classes used to store grids and data intended for visualization and/or checkpoint/restart include the term `image,' which evokes both colloquial and computing connotations (respectively `picture' and `a computer file capturing the contents of some computing element---hard drive, system memory, etc.').

\paragraph{FileSystemBasics}

In addition to some lower-level functionality, the classes in this division include two kinds of streams. 

{\tt ParametersStreamForm} is used to read scalar or one-dimensional array parameters of various datatypes from a text file, each line of which contains an entry of the form {\tt [name]=[value]} or {\tt [name]=[value1],[value2],...}. 
Units may be attached to real values with a tilde (\textasciitilde) followed by an all-caps string denoting one of the units available in {\tt UNIT\_Singleton} (see ``Specifiers'' in Section \ref{subsubsec:VariableManagement}).
For example, the default dimensionless value $1.0$ of the {\tt Amplitude} member of  
{\tt SineWaveDerivativeForm} (see Listings~\ref{lst:Module_SWD_Outline} and \ref{lst:Subroutine_Initialize_SWD}) could be overwritten with an entry like
\begin{lstlisting}
Amplitude=5.0~TESLA
\end{lstlisting} 
in a parameter file.
Parameters need not appear in the file in the order searched for by the program; indeed the absence of a parameter, or even the absence of the parameter file, will not in themselves cause a program exit or crash (though the lack of sensible default values for absent parameters certainly might).  
Our example program {\tt SineWaveDerivative} does not directly invoke objects of {\tt ParametersStreamForm}, but instead uses a higher-level parameter search facility described in Section \ref{subsubsec:Runtime} that looks both in a parameter file and at command line options. 

A second type of stream in {\tt FileSystemBasics} is {\tt TableStreamForm}, which reads simple two-dimensional tables of data from a text file.

\paragraph{GridImageBasics}

The remaining divisions in {\tt FileSystem}, beginning with {\tt GridImageBasics}, provide a fa\c{c}ade facilitating interaction with sophisticated I/O libraries. 
Different libraries may be implemented as `back ends', but the only specific implementation at present is the Silo library,\footnote{https://wci.llnl.gov/codes/silo} which is convenient for use with the visualization package VisIt.\footnote{https://wci.llnl.gov/codes/visit}  

The major class of interest in this division is {\tt GridImageStreamForm}, which handles some basic aspects of I/O that are the same regardless of the type(s) of grid(s) to be written. 
Its members and methods allow one to, for instance, open (for reading or writing) or close a file; make a directory (virtual within a file or actual in the file system, which may depend on the I/O library back end); list an existing file's contents; and automatically track filename numbering for consecutive write operations (e.g. for periodic output in time-dependent evolution).
In our example program, the method {\tt WriteVariables} of {\tt WaveDerivativeTemplate} declares and initializes an instance of {\tt GridImageStreamForm}---and uses that object to open a file---as follows:  
\begin{lstlisting}
  type ( GridImageStreamForm ) :: GIS

  call GIS
  call GIS
\end{lstlisting}
Here {\tt Name} and {\tt C} are aliases for the program's name and its global communicator. 
The argument used in calling the {\tt Open} method specifies the mode of access to the stream object. 
In this example, we ask the stream object to create a new file. 
With this access mode, a file number counter associated with the stream object is incremented and appended to the {\tt Name} argument to form the filename, and then a new file is created with this filename. 
This access mode also implies read and write access to the newly created file.
Other possible access modes are {\tt ACCESS\_WRITE}, which opens an existing file for writing by appending new data to the end of the file; and {\tt ACCESS\_READ}, which opens an existing file for reading only.\footnote{These access modes are parameter class members of the stream object (see also \footnotemark[9]).}
After the data is written, the call
\begin{lstlisting}
  call GIS
\end{lstlisting}
flushes the data associated with the stream\footnote{Depending on the I/O library and hardware, a call to the write method may only buffer the data and not actually write to the file system until the stream is closed.} and closes the file.

\paragraph{CurveImages} The class {\tt CurveImageForm} is used write a one-dimensional grid and data on it, which can be used to generate curves in a Cartesian $xy$ plane.
In our sample program the wave and its derivative (which can be used as dependent variables in curve plots) and the spatial coordinates of the lattice points (which can be used as the independent variable) are written with the following lines, which (aside from the declaration of the {\tt CurveImageForm} object {\tt CI}) lie between the calls to {\tt GIS \% Open} and {\tt GIS \% Close} just mentioned above:
\begin{lstlisting}
  type ( CurveImageForm ) :: &
    CI

  call CI
  call CI
         ( Directory = 'Curves', NodeCoordinate = Coordinate, nProperCells = WD
           oValue = 0, CoordinateUnitOption = CoordinateUnit )
  call CI
  call CI
\end{lstlisting}
The {\tt GridImageStream} object {\tt GIS} mentioned above is needed as an argument to the {\tt CI \% Initialize} method. 
The call to {\tt CI \% SetGrid} specifies a directory, gives the spatial coordinates of the lattice points, and includes an optional argument specifying the units of those coordinates.
The variable group {\tt VG\_Curves}, which comprises only the data of the wave and its derivative, was discussed in the ``VariableGroups'' portion of Section \ref{subsubsec:VariableManagement}.
Its use as an argument to the {\tt CI \% AddVariableGroup} method means that the variable data and metadata (names, units, etc.) specified by the variable group are available to the {\tt CurveImageForm} methods and can be included in the file as supported by the I/O library back end, which in turn allows a visualization package that supports the format to display it directly.
Here we attach only a single variable group, but any number can be added.
With these preparations, all that is needed is the call to {\tt CI \% Write}. 
Plots of the saved data---in which the previously specified units and variable names are automatically used by the visualization package VisIt without further intervention by the user---are shown in Figure~\ref{fig:SineWaveDerivative}.

\paragraph{StructuredGridImages, UnstructuredGridImages} 

The classes {\tt StructuredGridImageForm} and {\tt UnstructuredGridImageForm} are used to write two- or three-dimensional grids and data on them, which can be used to generate various kinds of multidimensional plots.
A structured grid may be rectilinear or curvilinear, but is in either case `logically rectangular': 
the quadrilateral or cuboid cells filling the domain can be indexed with a pair or triplet of integers, with a single list of coordinates in each dimension sufficing (together with the integer values of an index pair or triplet) to specify the edge values of any particular cell.
In an unstructured grid---which can be more irregular in shape---the coordinates of each individual node must be given, together with connectivity information that specifies how the nodes are arranged into cells.
The simple example problem in this section does not have any examples of structured or unstructured grids. 
Conceptually, however, their use follows that exemplified above in connection with curves: 
a {\tt SetGrid} method provides the information needed to set up a particular kind of grid, and an {\tt AddVariableGroup} method attaches dependent variables to the grid.
As we shall see in Section \ref{sec:Manifolds}, in the multilevel grids we use to represent charts (coordinate patches) in a manifold, the coarsest level is a structured grid and the refined levels are unstructured grids. 


\subsubsection{Runtime}
\label{subsubsec:Runtime}

The last division of {\tt VariableManagement} is the leaf division {\tt Runtime} (see the diagram on the left side of Figure \ref{fig:Basics_Structure}), which provides some generic functionality associated with running programs.
Included are a function that returns wall time, a command that displays memory usage, and a class that handles command line options (newly accessible to programs in Fortran 2003).

The user can access all of this functionality through the object {\tt PROGRAM\_HEADER} declared in {\tt PROGRAM\_HEADER\_Singleton}, which requires initialization.
Our example program {\tt SineWaveDerivative} contains the lines
\begin{lstlisting}
  allocate ( PROGRAM_HEADER )
  call PROGRAM_HEADER
         ( 'SineWaveDerivative', AppendDimensionalityOption = .false. )
\end{lstlisting}
which appear at line~9 of Listing~\ref{lst:Program_SWD_Outline}.
The first argument to the {\tt Initialize} method is the program name.
Because this program is strictly one-dimensional, a second optional argument overrides the default behavior of appending a short string indicating the dimensionality (e.g. {\tt '\_1D'}, {\tt '\_2D'}, etc.) to the program name.
Because most applications of GenASiS\ are expected to be written rather generically for anywhere from one to three space position space dimensions, the default behavior is to determine the dimensionality from a parameter file or command line option, or failing to find these, default to three dimensions. 
{\tt PROGRAM\_HEADER} has the global (`world') communicator as one of its members, and intialization of this singleton takes care of initialization of MPI (see Section~\ref{subsubsec:MessagePassing}).
 The {\tt UNIT} (see ``Specifiers'' in Section~\ref{subsubsec:VariableManagement}) and {\tt CONSOLE} (see Section~\ref{subsubsec:Display}) singletons are also initialized.
Once {\tt PROGRAM\_HEADER} is intialized its method {\tt ShowStatistics} can be called at any time to display total elapsed time, the amount of time spent in I/O, and current memory usage.

Another {\tt PROGRAM\_HEADER} method of significant interest is {\tt GetParameter}.
Our general approach to parameters\footnote{In this context we use the term `parameter' to mean a user-specified value that characterizes a problem or specifies some aspect of program execution. 
This usage is not quite as restrictive as the precise meaning of the Fortran {\tt parameter} keyword; but most of the time it ends up obeying the same spirit in practice, i.e. denoting some value that remains unchanged during program execution.} is to (1) set a reasonable default value in the code itself, (2) look for an overriding user-specified value in a parameter file, and (3) look for a higher-priority user-specified value among any command line options.
The {\tt PROGRAM\_HEADER \% GetParameter} method follows this philosophy:
it informs the user of the existing (presumably default) value, and then sequentially looks for and resets that value if it is found in a parameter file and/or on the command line (with appropriate notifications along the way).
The parameter file searched by default is one whose filename consists of the program name with the suffix {\tt '\_Program\_Parameters'} appended; 
alternatively, a {\tt ParametersStream} object (see ``FileSystemBasics'' in Section~\ref{subsubsec:FileSystem}) associated with some other parameter file can be supplied as an optional argument.
Our example program {\tt SineWaveDerivative} contains the line
\begin{lstlisting}
    call PROGRAM_HEADER
           ( SWD
\end{lstlisting}
at line~8 of Listing~\ref{lst:Subroutine_Initialize_SWD}.
The final optional argument returns the units that may (optionally) be specified by the parameter input (for instance in an expression like {\tt Amplitude=5.0{\textasciitilde}TESLA}). 


\section{MANIFOLDS}
\label{sec:Manifolds}

{\tt Manifolds} is the only division of {\tt Mathematics} that we report on in this paper.
The {\tt Mathematics} division of GenASiS\ contains generic mathematical constructs and solvers that are as agnostic as possible with regard to the specifics of any particular system, in order that they might in principle be used for a wide variety of physical (or other) problems.
(For the place of {\tt Mathematics} in the overall scheme of GenASiS, see Figure~\ref{fig:GenASiS_Structure}.) 
As illustrated in the left diagram of Figure~\ref{fig:Mathematics_Structure}, {\tt Manifolds} is subdivided into {\tt Charts} and {\tt Atlases}.
The middle diagram of Figure~\ref{fig:Mathematics_Structure} shows that {\tt Charts} includes the divisions {\tt ChartBasics}, {\tt Cells}, {\tt Submeshes}, {\tt Meshes}, {\tt Calculus}, and {\tt Intermeshes}.
All of these are leaf divisions except {\tt Calculus},
which includes the leaf divisions {\tt Differences} and {\tt Gradients} as shown in the right diagram of Figure \ref{fig:Mathematics_Structure}.

As was done in the case of {\tt Basics} in Section \ref{sec:Basics}, we illustrate some of the functionality available through {\tt Manifolds}---and continue to demonstrate our use of the object-oriented approach---in the context of an example problem.

\subsection{Another Problem and an Object-Oriented Solution}
\label{subsec:AnotherProblem}

The problem we use to illustrate the functionality of {\tt Manifolds} is somewhat similar to that discussed in Section \ref{sec:Basics}.
It includes the following tasks:
(a) initialize a manifold with a single `chart' or coordinate patch;
(b) refine the chart into a series of concentric spheres with increasing resolution;
(c) initialize a Gaussian function; 
(d) compute the gradient of this function;
(e) write the output in a format suitable for visualization.
We have contrived a Fortran {\tt program} (Listing~\ref{lst:Program_GG}) and two supporting {\tt module}s (Listings~\ref{lst:Module_GG_Outline} and \ref{lst:Module_FG_Outline}) to accomplish these tasks using the GenASiS\ {\tt Manifolds} and {\tt Basics} classes. 
\begin{lstlisting}[float,frame=tb,numbers=left,numbersep=5pt,xleftmargin=10pt,label=lst:Program_GG,caption={{\tt program GaussianGradient}}]
program GaussianGradient

  use Basics
  use GaussianGradient_Form  !-- See Listing \ref{lst:Module_GG_Outline}

  implicit none

  type ( GaussianGradientForm ) :: GG  !-- See Listing \ref{lst:Module_GG_Outline}

  allocate ( PROGRAM_HEADER )
  call PROGRAM_HEADER

  call GG
  call GG
  call GG
  call GG
  call GG

  deallocate ( PROGRAM_HEADER )

end program GaussianGradient
\end{lstlisting}

\begin{lstlisting}[float,frame=tb,numbers=left,numbersep=5pt,xleftmargin=10pt,label=lst:Module_GG_Outline,caption={Outline of {\tt module GaussianGradient\_Form}. Used at line~4 of Listing~\ref{lst:Program_GG}.}]
module GaussianGradient_Form

  use Basics
  use FunctionGradient_Template
  implicit none
  private

  type, public, extends ( FunctionGradientTemplate ) :: GaussianGradientForm
    type ( MeasuredValueForm ) :: Amplitude
    type ( MeasuredValueForm ), dimension ( 3 ) :: HalfWidth
  contains
    procedure, private, pass :: Initialize_GG
    generic :: Initialize => Initialize_GG
    procedure, private, pass :: Function
  end type GaussianGradientForm

contains

  !-- Definitions of {\tt subroutine Initialize\_GG} and {\tt function Function} omitted

end module GaussianGradient_Form

\end{lstlisting}

\begin{lstlisting}[float,frame=tb,numbers=left,numbersep=5pt,xleftmargin=10pt,label=lst:Module_FG_Outline,caption={Outline of {\tt module FunctionGradient\_Template}. Used at line~4 of Listing~\ref{lst:Module_GG_Outline}.}]
module FunctionGradient_Template

  use Basics
  use Manifolds
  implicit none
  private

  type, public, abstract :: FunctionGradientTemplate
    !-- Declaration of members omitted
  contains
    procedure, private, pass :: Initialize_FG       !-- See Listing~\ref{lst:SubroutineInitialize_FG}
    generic, public :: Initialize => Initialize_FG
    procedure, public, pass :: CreateNestedSpheres  !-- See Listing~\ref{lst:SubroutineCreateNestedSpheres}
    procedure, public, pass :: SetFunction          !-- See Listing~\ref{lst:SubroutineSetFunction}
    procedure, public, pass :: ComputeGradient      !-- See Listing~\ref{lst:SubroutineComputeGradient}
    procedure, public, pass :: Write                !-- See Listing~\ref{lst:SubroutineWrite}
    procedure ( FunctionInterface ), private, pass, deferred :: Function
  end type FunctionGradientTemplate

  !-- Declaration of the {\tt interface} for {\tt function FunctionInterface} omitted

contains

  !-- Definitions of {\tt subroutine}s {\tt Initialize\_FG}, {\tt CreateNestedSpheres},  {\tt SetFunction}, 
  ! \  \   {\tt ComputeGradient}, and {\tt Write} omitted

end module FunctionGradient_Template
\end{lstlisting}

The high-level structure of this example program is very similar to that of Section \ref{sec:Basics}.
The drivers {\tt SineWaveDerivative} (Listing~\ref{lst:Program_SWD_Outline}) and {\tt GaussianGradient} (Listing~\ref{lst:Program_GG}) each declare a single object and call its methods to carry out the required tasks. 
These objects each specialize a rather generic class to a particular functional form: {\tt SineWaveDerivative\_Form} (Listing~\ref{lst:Module_SWD_Outline}) extends {\tt WaveDerivative\_Template} (Listing~\ref{lst:Module_WD_Outline}), and {\tt GaussianGradient\_Form} (Listing~\ref{lst:Module_GG_Outline}) extends {\tt FunctionGradient\_Template} (Listing~\ref{lst:Module_FG_Outline}).
Beyond the members inherited from the parent class, these extensions add members parametrizing a specific functional form. 
The child classes also overload the {\tt Initialize} method and override the (deferred) {\tt Waveform} or {\tt Function} methods that return values of the specific functional form.
With the exception of {\tt Initialize}, the methods called by the driver {\tt program}s are fully declared and defined in the parent classes ({\tt WaveDerivative\_Template} and {\tt FunctionGradient\_Template}) and are agnostic with regard to the particular functional form whose derivative or gradient is to be taken.
These procedures can be called in the drivers as methods of the specialized extensions ({\tt SineDerivative\_Form} and {\tt GaussianGradient\_Form}) because they are inherited from the parent class.

There are also some differences between these example programs.
That the {\tt PROGRAM\_HEADER} initialization in line~12 of Listing~\ref{lst:Program_GG} for the {\tt GaussianGradient} program contains only one argument (the program name), while that for the {\tt SineWaveDerivative} program has two arguments (see Section~\ref{subsubsec:Runtime}), is a thread on which we can pull to unravel the distinctions.
Not having the optional argument {\tt AppendDimensionalityOption} present and set to {\tt .false.}, a command line option allows the {\tt GaussianGradient} program to be run in one dimension ({\tt Dimensionality=1D}), two dimensions ({\tt Dimensionality=2D}), or---and by default if no command line option is present---three dimensions ({\tt Dimensionality=3D}).
In contrast, the {\tt SineWaveDerivative} program is hard coded for one dimension only.
The {\tt GaussianGradient} program can accomplish its tasks not only in anywhere from one to three dimensions, but also with mesh refinement.
That it can do so with a line count comparable to that of the {\tt SineWaveDerivative} program is only possible because it draws on functionality available through {\tt Manifolds} classes: 
tasks such as setting up a grid, performing ghost cell data exchanges, and computing derivatives are easily done in {\tt SineWaveDerivative} with {\tt Basics} classes when only a single-level one-dimensional grid is involved, but {\tt Manifolds} provides higher-level functionality for performing such tasks in multidimensional and multilevel contexts with comparable ease.
A final difference we note is that {\tt SineWaveDerivative} uses a periodic lattice while {\tt GaussianGradient} has a boundary. 
Periodic boundary conditions are the default when using {\tt Manifolds} classes, but a more centralized function with a fixed boundary condition better suits our purpose of demonstrating a multilevel grid of a sort suitable for gravitational collapse. 


\subsection{Illustrations of {\tt Manifolds} Functionality}
\label{subsec:ManifoldsFunctionality}

With this example program in mind we describe some of the functionality available in {\tt Manifolds}. In particular we discuss in bottom-up fashion (i.e. in order of compilation) the code divisions appearing in Figure \ref{fig:Mathematics_Structure}.

\subsubsection{Charts}
\label{subsubsec:Charts}

{\tt Charts} is the first code division under {\tt Manifolds} in the diagram on the left side of Figure \ref{fig:Mathematics_Structure}.  
The term `chart' is mathematical jargon for a coordinate patch that covers an individual region of a manifold. 
We approximate the ideal of continuity with a finite sequence of meshes which provide, as necessary, increasing refinements of the coarsest (top-level) mesh. 
Our charts can be  one-, two-, or three-dimensional. 
Its subdivisions {\tt ChartBasics}, {\tt Cells}, {\tt Submeshes}, {\tt Meshes}, {\tt Calculus}, and {\tt Intermeshes} are shown in the middle diagram of Figure~\ref{fig:Mathematics_Structure}. 
We discuss each in turn.

\paragraph{ChartBasics}

The first division of {\tt Charts} is {\tt ChartBasics}, which contains some basic definitions needed by the subsequent classes that culminate in the class {\tt ChartForm} (see ``Intermeshes'' below). 
A number of parameters are defined in singleton objects; a few examples will appear later.

We note two classes in {\tt ChartBasics} that are used in connection with task (b) of our example program: refine the chart into a series of concentric spheres with increasing resolution.
{\tt AdaptionForm} is an extension of {\tt VariableGroupForm} (see ``VariableGroups'' in Section \ref{subsubsec:VariableManagement}) that contains variables related to mesh refinement and coarsening. 
{\tt GeometryBasicForm} is a different extension of {\tt VariableGroupForm} that contains the most basic geometric variables associated with the cells of a mesh (edge and center positions, face areas, volume, etc.) 
Instances of both these classes are referenced in a {\tt subroutine} called {\tt  SetRefineFlagLevel}, which marks a region for refinement; see Listing~\ref{lst:SubroutineSetRefineFlagLevel}.
(As we will see below, this routine for a single level is called at line~24 of Listing~\ref{lst:SubroutineCreateNestedSpheres} for the method {\tt CreateNestedSpheres} of {\tt FunctionGradientTemplate}.)
\begin{lstlisting}[float,frame=tb,numbers=left,numbersep=5pt,xleftmargin=10pt,label=lst:SubroutineSetRefineFlagLevel,caption={Outline of {\tt subroutine SetRefineFlagLevel}. Called at line~24 of Listing~\ref{lst:SubroutineCreateNestedSpheres}.}]
subroutine SetRefineFlagLevel ( FG, RefinementRadius, iLevel )

  class ( FunctionGradientTemplate ), intent ( inout ) :: FG
  real ( KDR ), intent ( in ) :: RefinementRadius
  integer ( KDI ), intent ( in ) :: iLevel

  integer ( KDI ) :: iV_1, iV_2  !-- iValues

  !-- Statements omitted
    
  associate &
    ( X => G
      Y => G
      Z => G
        
  associate ( R => sqrt ( X ** 2  +  Y ** 2  +  Z ** 2 ) )
  where ( R <= RefinementRadius )
    A
  end where
  end associate  !-- R

  end associate  !-- X, Y, Z
  !-- Other end associate statements omitted  
 
  end subroutine SetRefineFlagLevel
\end{lstlisting}
The argument {\tt iLevel} (lines~1 and 5) specifies the level of the multilevel grid structure on which this routine will act.
Lines~11-14 specify rank-one array aliases {\tt X}, {\tt Y}, and {\tt Z} for the coordinates of cell centers. 
{\tt G} is an alias to the instance of {\tt GeometryBasicForm} for this level. 
As an extension of {\tt VariableGroupForm}, it has a data array member {\tt Value}; the indices {\tt iV\_1} and {\tt iV\_2} specify a range of cells, and the elements of the array member {\tt G \% CENTER} index the columns of the data array containing the cell center values.
Lines~16-20 specify the spherical region targeted for refinement.
{\tt R} is an alias for the distance from the cell centers to the origin. 
{\tt A} is an alias to the instance of {\tt AdaptionForm} for this level.
Those cells for which {\tt R} is less than the argument {\tt RefinementRadius} (lines~1 and 4) have the `refine flag' variable---i.e. the column of {\tt A \% Value} specified by the member {\tt A \% REFINE\_FLAG}---set to a value (i.e. 1.0) that will trigger refinement. 

The last class in {\tt ChartBasics} that we mention is {\tt ChartHeaderForm}, which contains some basic parameters and other data describing a chart. 
Examples include the number of dimensions, the number of cells in each dimension on the coarsest level, the minimum and maximum coordinates in each dimension on the coarsest level, the number of levels, and so on.
The class {\tt ChartForm} (see ``Intermeshes'' below) is an extension of {\tt ChartHeaderForm}, so that all the members of {\tt ChartHeaderForm} are also members of {\tt ChartForm}.

\paragraph{Cells}

The second division of  {\tt Charts} is {\tt Cells}, which contains classes that enable construction of the adaptive structure that underlies our approximate representation of a three-dimensional continuous chart. 
This is an oct-tree (or, in restricted use in two dimensions or even one dimension, a quad- or binary tree respectively) that enables cell-by-cell refinement. 
The fundamental unit of the structure is a single finite cell, which is a segment, quadrilateral, or cuboid that can be split into two, four, or eight cells in one, two, or three dimensions respectively. 
A cell is represented by the class {\tt CellForm}, which contains pointers to its parent cell and to arrays of children and sibling cells. 
{\tt Cells} also contains the first of several interfaces to this oct-tree, a class {\tt CellListForm} that allows the creation of `cell lists,' or linked lists of cells. This can be used to, for example, create lists of selected parent cells in order to facilitate loops over frequently addressed subsets of cells on a particular level of the tree. 
Our example program does not directly reference any classes in {\tt Cells}.

\paragraph{Submeshes}

The next division of  {\tt Charts} is {\tt Submeshes}, which contains classes that define the next layer of interface to the oct-tree, a `submesh.' 
This entity, characterized by the class {\tt SubmeshForm}, is a grid in its own right; it consists of a subset of cells at a single level of the oct-tree, whose combined arrangement may be irregular in shape and even consist of multiple disconnected pieces. 
Among the members of {\tt SubmeshForm} are two cell lists, for `proper' and `ghost' cells, which allow the submesh to be domain-decomposed for parallel processing in a distributed-memory environment.
(Proper cells are typically the normal working computational cells assigned to a particular process. Ghost cells typically compose a partial boundary layer around the proper cells; 
their data is obtained through message passing with the neighboring processes that own those cells. 
See for instance Figure~\ref{fig:Decomposition}.) 
There are also members {\tt nProperCells} and {\tt nGhostCells} that store the number of proper and ghost cells, as well as a member {\tt oCell} giving the offset in a data array where storage begins for the subset of cells at a particular level composing the submesh.
These members are used, for example, to set the range of cells in Listing~\ref{lst:SubroutineSetRefineFlagLevel} discussed in ``ChartBasics'' above. 
In particular, the following lines are two of the statements omitted at line~9:
\begin{lstlisting}
  iV_1 = S
  iV_2 = S
\end{lstlisting}
Here {\tt S} is an alias to the instance of {\tt SubmeshForm} comprising the body of cells interior to the boundary---either the  boundary of the chart as a whole, or a coarse/fine boundary---at a particular level of refinement.
With these assignments, the expression {\tt iV\_1 : iV\_2} in lines~12-14 and 18 spans the range of data corresponding to all the proper and ghost cells in this submesh.
A submesh and data defined on it can independently be written to the file system, most often as an unstructured mesh due to its potentially irregular shape and even multiple pieces.
Examples of spherical meshes in two and three dimensions are shown in Figure~\ref{fig:Decomposition}.


\paragraph{Meshes}

The classes in the {\tt Meshes} division of {\tt Charts} combine multiple submeshes into a `mesh,' characterized by the class {\tt MeshForm}. 
A mesh includes four submeshes.
This is illustrated in Figure~\ref{fig:Submeshes}; shown in both two and three dimensional versions are the submeshes of the Level 1 (coarsest) mesh, which is rectangular, and the submeshes of the Level 2 (first refinement) mesh, which has a spherical shape. 
Two of these submeshes together comprise all the cells at a particular level of the oct-tree: the `Interior' submesh (or simply `Interior') contains all the normal computational cells, and the `Exterior' submesh includes all the cells that form a boundary layer---either the  boundary of the chart as a whole, or a coarse/fine boundary at the edge of a particular level---surrounding the Interior. 
The other two are the `Parents' and `Children' submeshes, which link the cells of the Interior to the adjacent coarser and finer levels respectively. 
The reason for the redundancies of a Level $i$ Parents with a Level $i-1$ Interior, and a Level $i$ Children with a Level $i+1$ Interior, is that the Interior submeshes on each level are independently domain-decomposed: the Level $i$ Parents and Children follow the decomposition of the Level $i$ Interior with which they are associated, which in general may be different from the decompositions of the Interior submeshes at Levels $i-1$ and $i+1$. 

A mesh is the typical locus of data storage. The storage available in an instance of {\tt VariableGroupForm} or its extensions (see ``VariableGroups'' in Section \ref{subsubsec:VariableManagement}) is subdivided into sections reserved for data associated with the cells of the four submeshes. 
A given cell in the oct-tree can play any or all of three roles---`level' (as part of an Interior or Exterior submesh), `parent' (as part of a Parents submesh), and `child' (as part of a Children submesh)---and as appropriate is assigned a number for each of these roles that corresponds to a row number of the data storage array in an instance of {\tt VariableGroupForm} or its extensions.

Task (c) of our example program for this section---initialize a Gaussian function---provides a further opportunity to see how this mesh-based storage is used.
(We have already caught initial glimpses in Listing~\ref{lst:SubroutineSetRefineFlagLevel} and the discussions thereof in ``ChartBasics'' and ``Submeshes'' above.)

Consider the execution of task (c) for a single mesh at a given level of refinement in {\tt subroutine SetFunctionLevel} in Listing~\ref{lst:SubroutineSetFunctionLevel}. 
(As we will see below, this routine for a single level is called at line~12 of Listing~\ref{lst:SubroutineSetFunction} for the method {\tt SetFunction} of {\tt FunctionGradientTemplate}.)
\begin{lstlisting}[float,frame=tb,numbers=left,numbersep=5pt,xleftmargin=10pt,label=lst:SubroutineSetFunctionLevel,caption={{\tt subroutine SetFunctionLevel}. Called at line~12 of Listing~\ref{lst:SubroutineSetFunction}.}]
subroutine SetFunctionLevel ( FG, iLevel )

  class ( FunctionGradientTemplate ), intent ( inout ) :: FG
  integer ( KDI ), intent ( in ) :: iLevel

  !-- associate statement omitted
  
  call FG
         ( [ nMaxCells, 1 ], NameOption = 'Variable', &
           VariableOption = [ FG
           UnitOption = [ FG

  call SetFunctionSubmesh ( FG, iLevel, SUBMESH
  if ( iLevel == 1 ) &
    call SetFunctionSubmesh ( FG, iLevel, SUBMESH

  !-- end associate statement omitted
  
end subroutine SetFunctionLevel
\end{lstlisting}
One of the members of {\tt FunctionGradientTemplate} declared at line~9 of Listing~\ref{lst:Module_FG_Outline} is an array of {\tt VariableGroupForm}, each element of which stores, for a given level of refinement, the variable whose gradient is to be taken: 
\begin{lstlisting}
  type ( VariableGroupForm ), dimension ( : ), allocatable :: Variable
\end{lstlisting}
In Listing~\ref{lst:SubroutineSetFunctionLevel}, the argument {\tt iLevel} (lines~1 and 4) specifies which level is to have the variable set.
The element of {\tt FG \% Variable} for this level is initialized in lines~8-11; for more on initialization of instances of {\tt VariableGroupForm}, see ``VariableGroups'' in Section \ref{subsubsec:VariableManagement}. (The value {\tt nMaxCells} is an alias to the maximum number of cell data values instance of {\tt MeshForm} can have.)
Separate calls to another {\tt subroutine}, {\tt SetFunctionSubmesh} shown in Listing~\ref{lst:SubroutineSetFunctionSubmesh}, set the variable whose gradient is to be taken on the Interior submesh (line~13) and, for the coarsest level, the Exterior submesh as well (lines~14-15). Values on the Exterior of the coarsest level are set here because this constitutes the boundary of the chart as a whole; see the left panels of Figure~\ref{fig:Submeshes}. 
In a physics problem in which a single chart covered the entire computational domain, the Exterior submesh on the coarsest level is where physical boundary conditions would be applied; here we just apply a specified functional form whose gradient is to be taken.
For the multilevel grid configuration in our example problem in this section, the Exterior submeshes of finer levels are not on the boundary on the chart as a whole; see for instance the right panels of Figure~\ref{fig:Submeshes}. 
In our simple example problem, values for the Exterior on the finer levels could be directly specified by the known functional form, but to illustrate what would happen in a more realistic problem, we instead set these Exterior values on finer levels later by interpolation from the coarser level (i.e. by `prolongation'; see ``Intermeshes'' below).
Finally, we note that the parameters {\tt SUBMESH \% INTERIOR} and {\tt SUBMESH \% EXTERIOR} appearing in lines~13 and 15 are members a singleton {\tt SUBMESH}; this is one of the singleton objects mentioned in ``ChartBasics'' above.

We turn next to the execution of task (c) for a single submesh, shown in {\tt subroutine SetFunctionSubmesh} in Listing~\ref{lst:SubroutineSetFunctionSubmesh} and called at lines~11 and 13 of Listing~\ref{lst:SubroutineSetFunctionLevel}.
\begin{lstlisting}[float,frame=tb,numbers=left,numbersep=5pt,xleftmargin=10pt,label=lst:SubroutineSetFunctionSubmesh,caption={Outline of {\tt subroutine SetFunctionSubmesh}. Called at lines~11 and 13 of Listing~\ref{lst:SubroutineSetFunctionLevel}.}]
subroutine SetFunctionSubmesh ( FG, iLevel, iSubmesh )

  class ( FunctionGradientTemplate ), intent ( inout ) :: FG
  integer ( KDI ), intent ( in ) :: iLevel, iSubmesh

  integer ( KDI ) :: iV_1, iV_2  !-- iValues

  !-- associate statements omitted

  iV_1 = S
  iV_2 = S
         
  associate &
    ( X => G
      Y => G
      Z => G
      V => FG

  V = FG

  end associate  !-- X, Y, Z
  !-- Other end associate statements omitted

  end subroutine SetFunctionSubmesh
\end{lstlisting}
Lines~10 and 11 are similar, but not identical, to the assignments shown in ``Submeshes'' above in connection with Listing~\ref{lst:SubroutineSetRefineFlagLevel}.
Here {\tt S} is an alias to the instance of {\tt SubmeshForm} indicated by the arguments {\tt iLevel} and {\tt iSubmesh} (lines~1 and 4).
Unlike the case seen previously, the assignments in lines~10 and 11 cause the range {\tt iV\_1 : iV\_2} in lines~14-17 to span only the proper cells in this submesh, and not the ghost cells. 
Despite the fact that we have a known functional form in this example problem, we do not fill in values for ghost cells here.
Instead, to illustrate what is typically needed in a more realistic problem, we instead set these ghost cell values later by an `exchange' of ghost cell data with neighboring processes (using message passing; see ``Calculus'' below).
Aliases for rank-one arrays for the coordinates of the centers of the proper cells, and for the values in those cells of the variable whose gradient is to be taken, are set up in lines~13-17.
The variable whose gradient is to be taken is assigned via the function call in line~19. 
As far as {\tt FunctionGradientTemplate} is concerned, {\tt FG \% Function} could be anything; its definition is {\tt deferred} and only an interface is specified in lines~17 and 20 of Listing~\ref{lst:Module_FG_Outline}.
The functional form is specified to be a Gaussian in {\tt GaussianGradientForm} outlined in Listing~\ref{lst:Module_GG_Outline}. 
This class extends {\tt FunctionGradientTemplate} and overrides the previous {\tt deferred} declaration, as indicated in lines~14 and 19 of Listing~\ref{lst:Module_GG_Outline}.

Along with being the typical locus of data storage, {\tt MeshForm} also has members with connectivity information, that is, lists of cell numbers and their siblings' cell numbers, in order to facilitate certain operations on particular sets of selected cells without having to walk through the oct-tree, access sibling pointers, etc. 
The application of boundary conditions is one example of an operation that uses this sort of connectivity information. 
Another is the overlapping of work and communication in time-explicit solves requiring only near-neighbor information (e.g. in hydrodynamics): `exchange' cells (those whose data must be communicated to the ghost cells of other message-passing processes) can be updated, and non-blocking communication initiated, before `non-exchange' cells are updated. 
Connectivity information for these two categories of cells allows data to be loaded from their non-contiguous loci in variable group storage into contiguous storage in a packed variable group (see ``VariableGroups'' in Section \ref{subsubsec:VariableManagement}) for efficient processing by the solver. 


\paragraph{Calculus}
 
The differential structure available in a chart is realized in discrete approximation in the classes defined in the {\tt Calculus} division of {\tt Charts}. 
As shown in the right diagram of Figure~\ref{fig:Mathematics_Structure}, this division includes such operations as {\tt Differences} and {\tt Gradients}, which are implemented on a single level of refinement in connection with the class {\tt MeshForm} just discussed above.

Task (d) of our example program for this section---compute the gradient of a function, in this case a Gaussian---utilizes the class {\tt GradientMeshForm} in {\tt Gradients}, and also provides an example of overlapping work and communication of ghost cell data.
Note that {\tt ComputeGradient} is declared as a method of {\tt FunctionGradientTemplate} in line~15 of Listing~\ref{lst:Module_FG_Outline}; this {\tt subroutine} is shown in Listing~\ref{lst:SubroutineComputeGradient}.
\begin{lstlisting}[float,frame=tb,numbers=left,numbersep=5pt,xleftmargin=10pt,label=lst:SubroutineComputeGradient,caption={{\tt subroutine ComputeGradient}. Fits in lines~24-25 of Listing~\ref{lst:Module_FG_Outline}. Called at line~16 of Listing~\ref{lst:Program_GG}.}]
subroutine ComputeGradient ( FG )

  class ( FunctionGradientTemplate ), intent ( inout ) :: FG

  integer ( KDI ) :: iL  !-- iLevel

  !-- associate statement omitted

  allocate ( FG

  do iL = 1, C

    associate &
      ( M  => C
        V  => FG
        GM => FG

    call GM
    call M 
    call GM
    call M 
    call GM

    end associate  !-- M, V, GM

  end do

  !-- end associate statement omitted

end subroutine ComputeGradient
\end{lstlisting}
One of the members of {\tt FunctionGradientTemplate} declared at line~9 of Listing~\ref{lst:Module_FG_Outline} is an array of {\tt GradientMeshForm}, each element of which contains the gradient for a given level of refinement: 
\begin{lstlisting}
 type ( GradientMeshForm ), dimension ( : ), allocatable :: Gradient
\end{lstlisting}
This member is allocated at line~9 of Listing~\ref{lst:SubroutineComputeGradient}.
Here {\tt C} is an alias to the instance of {\tt ChartForm} used in this program, and its member {\tt nLevels} is one of those inherited from {\tt ChartHeaderForm} as mentioned in ``ChartBasics'' above.
There follows a loop over the levels of the chart in lines~11-26. 
Aliases {\tt M}, {\tt V}, and {\tt GM} to the instances of the {\tt MeshForm}, of the {\tt VariableGroupForm} whose gradient is to be taken, and of the {\tt GradientMeshForm} for a particular level of refinement are set in lines~14-16.
The instance of {\tt GradientMeshForm} is initialized in line~18; in addition to a name, it needs the instance of {\tt MeshForm} with which it is associated, and the instance of {\tt VariableGroupForm} whose gradient is to be taken.
`Exchange' cells---those on the boundary between domains assigned to different processes---need the data in their ghost cell neighbors before a gradient can be taken.
This `non-blocking' exchange of ghost cell data is initiated in line~19. 
Meanwhile the gradient values in `non-exchange' cells, which do not have any ghost cell neighbors, can be computed in line~20 while the ghost cell data exchange continues in the background. 
Once the ghost data exchange is completed (line~21), the gradient values in exchange cells can be computed (line~22).
(Connectivity lists for exchange and non-exchange cells were discussed at the end of ``Meshes'' above. 
{\tt CONNECTIVITY} is one of the singleton objects mentioned in ``ChartBasics'' above; here its members {\tt EXCHANGE} and {\tt NONEXCHANGE} are parameters identifying the corresponding connectivity lists.)

\paragraph{Intermeshes}

The classes in the {\tt Intermeshes} division of {\tt Charts} bring together the structures described above into a `chart,' or coordinate patch, characterized by the class {\tt ChartForm}.
As already suggested above, {\tt ChartForm} has among its members an array of {\tt MeshForm}, each element of which corresponds to a (potential) level of the refinable oct-tree. 
This member of {\tt ChartForm}, an array of {\tt MeshForm}, is named {\tt Level}. 
Initialization of an object of class {\tt ChartForm} automatically includes setting up the coarsest level of the chart in accordance with parameters mentioned above in ``ChartBasics.''

Each successive level of the chart is a refinement of the previous level, allowing the chart to approximate the ideal of continuity as needed.
An example, generated by our sample problem for this section, is shown in Figure~\ref{fig:Chart2D}.
It consists of 10 levels and illustrates the dynamic range in length scales accessed during the gravitational collapse of the core of a massive star.
The levels of this chart are generated by the method {\tt CreateNestedSpheres} of {\tt FunctionGradientTemplate}, declared in line~13 of Listing~\ref{lst:Module_FG_Outline}; this {\tt subroutine} is shown in Listing~\ref{lst:SubroutineCreateNestedSpheres}.
\begin{lstlisting}[float,frame=tb,numbers=left,numbersep=5pt,xleftmargin=10pt,label=lst:SubroutineCreateNestedSpheres,caption={{\tt subroutine CreateNestedSpheres}. Fits in lines~24-25 of Listing~\ref{lst:Module_FG_Outline}. Called at line~14 of Listing~\ref{lst:Program_GG}.}].
subroutine CreateNestedSpheres ( FG )

  class ( FunctionGradientTemplate ), intent ( inout ) :: FG

  integer ( KDI ) :: iL, nLevels
  real ( KDR ), dimension ( : ), allocatable :: MeshRadius

  !-- associate statement omitted

  nLevels = 1
  call PROGRAM_HEADER
         ( nLevels, 'nLevels', IgnorabilityOption = CONSOLE

  if ( nLevels == 1 ) return

  allocate ( MeshRadius ( 2 : nLevels ) )
  MeshRadius ( 2 ) &
    = 0.5_KDR * minval ( C
  do iL = 3, nLevels
    MeshRadius ( iL ) = 0.5_KDR * MeshRadius ( iL - 1 )
  end do

  do iL = 1, nLevels - 1
    call SetRefineFlagLevel ( FG, MeshRadius ( iL + 1 ), iL )  !-- See Listing~\ref{lst:SubroutineSetRefineFlagLevel}
    call C
    call C
  end do

  !--  end associate statement omitted

end subroutine CreateNestedSpheres
\end{lstlisting}
The local variable {\tt nLevels} is given a default value of 1 (lines~5, 10). 
A larger value of {\tt nLevels} can be specified in a parameter file and/or the command line using the {\tt GetParameter} method of the {\tt PROGRAM\_HEADER} singleton (line~11; see Section~\ref{subsubsec:Runtime}).
If {\tt nLevels > 1} (line~14), the radii of spherical meshes at levels {2 $\le$ {\tt iL} $\le$ {\tt nLevels}} are computed (lines~16-21), where the local variable {\tt iL} (line~5) indexes the levels. 
The refinements that set up levels 2 through {\tt nLevels} occur in lines~23-27.
The routine {\tt SetRefineFlagLevel} called in line~24 was given in Listing~\ref{lst:SubroutineSetRefineFlagLevel} and discussed in ``ChartBasics'' earlier in this section.
Invocations of the methods {\tt AddLevel} and {\tt SetLevel} of {\tt ChartForm} (lines~25-26) prepare a new, refined level. 
{\tt C} is an alias to the instance of {\tt ChartForm} used in this program. 
The method {\tt AddLevel} performs the initializations and storage allocations needed for a new level, and increments the member {\tt C \% nLevels}.
The refinement actually takes place in the call to {\tt C \% SetLevel}, along with the setting of the geometry, boundary, and connectivity of the new level.

The interpolation of Level~$i$ data to Level~$i+1$ is called `prolongation,' and the averaging of Level~$i$ data to Level~$i-1$ is called `restriction.' 
These operations are illustrated in Figure~\ref{fig:ProlongationRestriction}.
The first step of prolongation (left panel in Figure~\ref{fig:ProlongationRestriction}) is to interpolate data from the Level~$i$ Interior submesh to the Level~$i$ Children submesh.\footnote{Gradients used in the interpolation optionally can be computed with a slope limiter in order to respect discontinuities.} 
The Level~$i$ Children submesh is refined relative to the Level~$i$ Interior with which it is associated.
But in a distributed-memory message passing environment, the Level~$i$ Children follows the domain decomposition of the Level~$i$ Interior, as mentioned in ``Meshes'' above, making the interpolation step a local operation. 
This is why both the Interior and Children submeshes are part of the same Level~$i$ mesh.
But the distribution of the Level~$i+1$ Interior cells among processes in general may differ from that of the Level~$i$ Interior, as the Interior submeshes on each level are independently domain-decomposed. 
Thus the second and final step of prolongation is communication, in general involving message passing, between from the Level~$i$ Children to the Level~$i+1$ Interior and/or Exterior. 
Similarly, restriction (right panel in Figure~\ref{fig:ProlongationRestriction}) involves first an averaging of Level~$i$ Interior and/or Exterior data to the Level~$i$ Parents (a local operation), followed by communication from the Level~$i$ Parents to the Level~$i-1$ Interior (which in general involves message passing).

In our example program, prolongation and restriction are performed in the {\tt SetFunction} method of {\tt FunctionGradientTemplate}.
This method is declared in line~14 of Listing~\ref{lst:Module_FG_Outline} and the {\tt subroutine} is shown in Listing~\ref{lst:SubroutineSetFunction}.
\begin{lstlisting}[float,frame=tb,numbers=left,numbersep=5pt,xleftmargin=10pt,firstnumber=1,label=lst:SubroutineSetFunction,caption={{\tt subroutine SetFunction}. Fits in lines~24-25 of Listing~\ref{lst:Module_FG_Outline}. Called at line~15 of Listing~\ref{lst:Program_GG}.}].
subroutine SetFunction ( FG )

  class ( FunctionGradientTemplate ), intent ( inout ) :: FG

  integer ( KDI ) :: iL  !-- iLevel

  !-- associate statement omitted

  allocate ( FG

  do iL = 1, C
    call SetFunctionLevel ( FG, iL )  !-- See Listing~\ref{lst:SubroutineSetFunctionLevel}
  end do

  do iL = C
    call C
           ( FG
             Interior = .true., Exterior = .false., iLevel = iL )
  end do

  do iL = 1, C
    call C
           ( FG
             Interior = .false., Exterior = .true., iLevel = iL )
  end do

  !-- end associate statement omitted

end subroutine SetFunction
\end{lstlisting}
The member {\tt Variable} of {\tt FunctionGradientTemplate}, allocated in line~9, 
is an array of {\tt VariableGroupForm} that we encountered in ``Meshes'' above.
Each element of the array stores, for a given level of refinement, the variable whose gradient is to be taken.
As before, {\tt C} is an alias to the instance of {\tt ChartForm} used in this program.
The routine {\tt SetFunctionLevel}, given in Listing~\ref{lst:SubroutineSetFunctionLevel} and discussed in ``Meshes'' above, is called for all levels in lines~11-13.
An example of restriction from the finest to the coarsest level follows in lines~15-19.
In this toy problem the example is gratuitous because the Interior at each level already has been initialized from a known analytic function. 
In a more realistic problem, however, restriction of more accurate data from a finer level to a coarser one can be a necessary operation.
In this particular case the arguments specify that only the Interior data on Level~$i$ is to be averaged to Level~$i-1$; see the right panel of Figure~\ref{fig:ProlongationRestriction}.
Finally, lines~21-25 illustrate prolongation from the coarsest to the finest level.
In this case data is interpolated from the Interior of Level~$i$ to the Exterior of Level~$i+1$, providing boundary data needed in this example to compute the gradient at the finer levels.
See the left panel of Figure~\ref{fig:ProlongationRestriction}.


\subsubsection{Atlases}
\label{subsubsec:Atlases}

{\tt Atlases} is the second and final code division under {\tt Manifolds} in the diagram on the left side of Figure \ref{fig:Mathematics_Structure}.  
The term `atlas' is mathematical parlance for a collection of charts that covers a manifold. 
A sufficiently simple manifold can be described with a single chart; this is the case with the {\tt GaussianGradient} example problem in this section. 
But one can imagine many reasons to use multiple charts. 
For example, a `yin-yang' atlas includes two overlapping charts that---like the two pieces of leather that form the surface of a baseball---cover a three-dimensional space with separate spherical coordinates in such a way as to avoid coordinate singularities at the origin and on the polar axis \citep[e.g.][]{Wongwathanarat2010b}. 
One could efficiently handle both spherically symmetric radial gravitational collapse and multdimensional phenomena at smaller radius by marrying a central three-dimensional Cartesian mesh to a one-dimensional radial mesh that begins near the surface of the Cartesian box and extends to much larger radius \citep{Scheidegger2008}. 
Choices of algorithm and approximation might suggest that different pieces of physics be treated on separate charts, with some facility for interpolation between them \citep{Rampp2002,Buras2006a}. 
A binary stellar system could be handled by covering the two bodies with separate meshes that move within a larger mesh to follow the orbital motion \citep{Scheel2006}. 
Some manifolds, such as some spaces described by general relativity, may be sufficiently complicated as to require multiple charts for a reasonable description. 
Another example is the phase space---position space plus momentum space---needed for relativistic kinetic theory, which in mathematical idealization contains an infinite number of charts: this is a `tangent bundle,' consisting of a (curved, in the general case) base space, i.e position space, together with a flat tangent space, i.e. a momentum space, at every point of the base space. 

The class {\tt AtlasForm} has among its members an array {\tt Chart} of {\tt ChartForm}.
In our example program, an instance of {\tt AtlasForm} is one of the members of {\tt FunctionGradientTemplate} declared at line~9 of Listing~\ref{lst:Module_FG_Outline}:
\begin{lstlisting}
  type ( AtlasForm ), allocatable :: Atlas
\end{lstlisting}
This member is initialized in the method {\tt Initialize\_FG} of {\tt FunctionGradientTemplate}, declared in line~11 of Listing~\ref{lst:Module_FG_Outline} and shown in Listing~\ref{lst:SubroutineInitialize_FG}.
\begin{lstlisting}[float,frame=tb,numbers=left,numbersep=5pt,xleftmargin=10pt,firstnumber=1,label=lst:SubroutineInitialize_FG,caption={{\tt subroutine Initialize\_FG}. Fits in lines~24-25 of Listing~\ref{lst:Module_FG_Outline}.}].
  subroutine Initialize_FG ( FG, FunctionName, FunctionUnit )

    class ( FunctionGradientTemplate ), intent ( inout ) :: FG
    character ( * ), intent ( in ) :: FunctionName
    type ( UnitForm ), intent ( in ) :: FunctionUnit

    FG
    FG

    allocate ( FG
    associate ( A => FG

    call A
    call A

    end associate  !-- A

  end subroutine Initialize_FG
\end{lstlisting}
Allocation of {\tt FG \% Atlas} occurs in line~10 and an alias {\tt A} to it is specified in line~11.
Initialization of the atlas is accomplished with the call in line~13.
After initialization one or more charts must be added.
This example involves only a single chart, which is added with a call to the method {\tt AddChart} in line~14.
The name of the chart, provided as an argument in this call, is used to find the parameter file needed to set up the chart.
In previous listings in this section we have used an alias {\tt C} to refer to the single chart used in our example program.
This alias is set with the expression
\begin{lstlisting}
  associate ( C => FG
\end{lstlisting}
which would appear for instance at line~7 of Listing~\ref{lst:SubroutineComputeGradient}, line~8 of Listing~\ref{lst:SubroutineCreateNestedSpheres}, and line~7 of Listing~\ref{lst:SubroutineSetFunction}.

{\tt AtlasForm} has methods for writing its charts, and data thereon, to output files.
These are called in the method {\tt Write} of {\tt FunctionGradientTemplate}, declared in line~16 of Listing~\ref{lst:Module_FG_Outline} and shown in Listing~\ref{lst:SubroutineWrite}.
\begin{lstlisting}[float,frame=tb,numbers=left,numbersep=5pt,xleftmargin=10pt,firstnumber=1,label=lst:SubroutineWrite,caption={{\tt subroutine Write}. Fits in lines~24-25 of Listing~\ref{lst:Module_FG_Outline}. Called at line~17 of Listing~\ref{lst:Program_GG}.}].
  subroutine Write ( FG )

    class ( FunctionGradientTemplate ), intent ( inout ) :: FG

    integer ( KDI ) :: iL, iD
    type ( GridImageStreamForm ) :: GIS

    associate ( A => FG

    allocate ( FG
    do iL = 1, C
      associate ( GM => FG
      do iD = 1, C
        call GO ( iD )
      end do
      end associate  !-- G, GO
    end do

    call GIS
    call GIS
    call A

    call A
    call A
    call A
    do iD = 1, C
      call A
    end do

    call A

    call A
    call GIS

    end associate  !-- A, C

  end subroutine Write
\end{lstlisting}
Aliases {\tt A} and {\tt C} for the atlas and its chart are set in line~8.
An instance of {\tt GridImageStreamForm} (lines~6, 19-20, 33; see ``GridImageBasics'' in Section~\ref{subsubsec:FileSystem}) provides access to an output file. 
The writing of the atlas to this output file is prepared by a call to the method {\tt OpenStream} of {\tt AtlasForm} (line~21), and completed by a call to the method {\tt CloseStream} (line~32). 
The fields (i.e. variables) to be written are specified by calls to the method {\tt AddFieldImage}.
Three of these, in lines~23-25, are straightforward: {\tt C \% Geometry}, {\tt C \% Adaption}, and {\tt FG \% Variable} (see ``ChartBasics'' and ``Meshes'' in Section~\ref{subsubsec:Charts}) are already arrays of {\tt VariableGroupForm} over the levels of the chart corresponding to the argument {\tt iChart = 1}.
Output of the gradient is a bit more complicated because the storage is not innately arranged as arrays over levels as expected by the {\tt AddFieldImage} method.
A rank-two array {\tt GradientOutput} of {\tt VariableGroupForm}, whose first index spans number of levels and whose second index spans the number of dimensions, is prepared in lines~10-17 and used in calls to {\tt AddFieldImage} in lines~26-28.
The write itself occurs with the call to the {\tt Write} method of {\tt AtlasForm} in line~30.

Plots made from data produced by our {\tt GaussianGradient} example program with {\tt nLevels=10} are presented in Figures~\ref{fig:Gaussian2D}-\ref{fig:GaussianGradient3D}. 
These show Gaussian functions in two and three dimensions and the $x$ component of their gradients, illustrating the range of scales encountered in the gravitational collapse of the core of a massive star and the continuity of the Gaussian function and its gradient at boundaries between adjacent levels of refinement. 


\section{CONCLUSION}
\label{sec:Conclusion}

This paper is the first in a series on GenASiS ({\em Gen}eral {\em A}strophysical {\em Si}mulation {\em S}ystem), a new code ultimately aimed at state-of-the-art simulations of core-collapse supernovae and other astrophysical problems.
Distinguishing features include an object-oriented approach and cell-by-cell mesh refinement.
The core of GenASiS includes three major divisions: 
{\tt Basics}, which contains some basic utilitarian functionality for large-scale simulations on distributed-memory supercomputers;
{\tt Mathematics}, which includes generic mathematical constructs and solvers that are as agnostic as possible with regard to the specifics of any particular system;
and {\tt Physics}, which sets up physical spaces associated with various theories of spacetime (including gravity), defines various forms of stress-energy, and combines these into `universes.'
As illustrated in Figure~\ref{fig:GenASiS_Structure}, these three divisions are present in each of the two highest-level divisions shown:
{\tt Modules}, which contains class implementations; and {\tt Programs}, which contains a unit test for each class.
{\tt Programs} also has a division {\tt Applications} containing drivers for applications beyond unit tests.
In this first paper we focus on {\tt Basics} and on the {\tt Manifolds} division of {\tt Mathematics}.
These contain fundamental capabilities that underpin implementations of the solvers and physics needed for the simulation of core-collapse supernovae and other astrophysical systems.

The content of {\tt Basics} is outlined in Figure~\ref{fig:Basics_Structure}.
The first division, {\tt VariableManagement}, contains several subdivisions shown in the right diagram of Figure~\ref{fig:Basics_Structure}.
{\tt Specifiers} contains classes used in the specification of number and character variables, as well as mathematical and physical constants and a means of dealing with units.
{\tt ArrayOperations} includes some basic operations on arrays.
{\tt ArrayArrays} contains classes that can be used to form arrays of arrays;
one application is the construction of so-called ragged arrays.
{\tt VariableGroups} contains classes we use extensively in handling collections of variables, especially sets of related physical fields. 
After {\tt VariableManagement} in the left diagram of Figure~\ref{fig:Basics_Structure} comes {\tt Display}, which contains infrastructure for displaying messages and variables to the standard output in a uniform and orderly way.
The classes in {\tt MessagePassing}, found in its subdivisions {\tt MessagePassingBasics}, {\tt PointToPoint}, and {\tt Collective} (see the middle right diagram in Figure~\ref{fig:Basics_Structure}), abstract the data and methods useful for working in a message passing parallel computing environment.
{\tt FileSystem} is next in the left diagram of Figure~\ref{fig:Basics_Structure}; 
the classes in {\tt FileSystemBasics}, {\tt GridImageBasics}, {\tt CurveImages}, {\tt StructuredGridImages}, and {\tt UnstructuredGridImages} (see the middle left diagram in Figure~\ref{fig:Basics_Structure}) handle I/O to disk, in part by providing a fa\c{c}ade facilitating interaction with sophisticated I/O libraries.
Finally, {\tt Runtime} provides a {\tt PROGRAM\_HEADER} object giving access to such facilities such as a function that returns wall time, a command that displays memory usage, and a class that handles command line options (newly accessible to programs in Fortran 2003).

The content of {\tt Manifolds}, a division of {\tt Mathematics}, is outlined in Figure~\ref{fig:Mathematics_Structure}.
The first division in the left diagram of Figure~\ref{fig:Mathematics_Structure}, {\tt Charts}, contains several subdivisions shown in the middle diagram of Figure~\ref{fig:Mathematics_Structure}; these comprise the infrastructure representing a coordinate patch that covers an individual region of a manifold. 
We approximate the ideal of continuity with a finite sequence of meshes which provide, as necessary, increasing refinements of the coarsest (top-level) mesh. 
Our charts can be  one-, two-, or three-dimensional. 
{\tt ChartBasics} contains some basic definitions needed by the subsequent classes that culminate in the class {\tt ChartForm}.
{\tt Cells} contains classes that enable construction of the adaptive structure that underlies our approximate representation of a three-dimensional continuous chart---an oct-tree (or, in restricted use in two dimensions or even one dimension, a quad- or binary tree respectively) that enables cell-by-cell refinement.
{\tt Submeshes} includes classes that culminate in {\tt SubmeshForm}, which embodies a grid consisting of a subset of cells at a single level of the oct-tree, whose combined arrangement may be irregular in shape and even consist of multiple disconnected pieces.
The class {\tt MeshForm} in {\tt Meshes} integrates four submeshes into a single level of the multilevel chart; 
the Interior and Exterior submeshes comprise all the cells at a particular level of the oct-tree, and the Parents and Children submeshes link the cells of the Interior to the adjacent coarser and finer levels respectively (see Figure~\ref{fig:Submeshes}). 
The differential structure available in a chart is realized in discrete approximation in the classes defined in {\tt Calculus}; as shown in the right diagram of Figure~\ref{fig:Mathematics_Structure}, this includes the operations {\tt Differences} and {\tt Gradients}, which are implemented on a single level of refinement in connection with the class {\tt MeshForm}. 
The classes in {\tt Intermeshes} culminate in {\tt ChartForm}, whose array of {\tt MeshForm} allows for a refinable chart or coordinate patch, with capabilities for interpolation of data to a more refined level (prolongation) and averaging of data to coarser level (restriction).
After {\tt Charts}, the second and final division in {\tt Manifolds} in the left diagram of Figure~\ref{fig:Mathematics_Structure} is {\tt Atlases}, which defines a collection of charts covering a manifold.

We present two sample problems, one each for {\tt Basics} and {\tt Manifolds}. 
The {\tt SineWaveDerivative} program, which illustrates the divisions of {\tt Basics},
computes the derivative of a sine wave in one dimension on a periodic lattice, with output shown in Figure~\ref{fig:SineWaveDerivative}.
The {\tt GaussianGradient} program, which illustrates the divisions of {\tt Manifolds},
computes the gradient of a Gaussian function in any of one to three dimensions on a multilevel grid of a sort suitable for gravitational collapse, with output shown in Figures~\ref{fig:Gaussian2D}-\ref{fig:GaussianGradient3D}.
That it can do so with a line count comparable to that of the {\tt SineWaveDerivative} program is only possible because it draws on functionality available through {\tt Manifolds} classes: 
tasks such as setting up a grid, performing ghost cell data exchanges, and computing derivatives are easily done in {\tt SineWaveDerivative} with {\tt Basics} classes when only a single-level one-dimensional grid is involved, but {\tt Manifolds} provides higher-level functionality for performing such tasks in multidimensional and multilevel contexts with comparable ease. 

Listings drawn from these sample programs serve not only to exhibit the contents of the {\tt Basics} and {\tt Manifolds} divisions of GenASiS, but also to illustrate our use of object-oriented design principles using the features of Fortran 2003 that support this programming paradigm \citep[e.g.][]{Reid2007,Adams2008}.
These principles promote code reusability by facilitating generality, extensibility, and maintainability.
In our context, an object-oriented approach enables the flexibility connoted by the `General' in GenASiS---the capacity of the code to include and refer to multiple algorithms, solvers, and physics and numerics choices with the same abstracted names and/or interfaces.
{\em Abstraction} identifies the major concepts required by a program, without specifying the details of implementation.
{\em Encapsulation} bundles together and controls access to the data (or {\em members}) and actions (or {\em methods}) associated with a particular concept into a self-contained unit---a {\em class}.
Done well, abstraction and encapsulation lead to {\em decoupling}---the separation of code into building blocks that are as independent and reusable as possible. 
Reusability is further enhanced by {\em polymorphism}, which enables multiple versions or implementations of the same basic concept to be referred to and used interchangeably.
Closely related is {\em inheritance,} which allows new ({\em child}) classes to be formed from prior ({\em parent}) classes: existing members and methods are retained or modified, and new members and methods can be added.

By way of applying the principles of abstraction, encapsulation, and decoupling, we have established conventions for implementing classes via Fortran {\tt module}s and derived {\tt type}s.
These conventions are illustrated in Listings~\ref{lst:Module_SWD_Outline} and \ref{lst:Module_WD_Outline} from the {\tt SineWaveDerivative} example, and Listings~\ref{lst:Module_GG_Outline} and \ref{lst:Module_FG_Outline} from the {\tt GaussianGradient} example.
A single {\tt module} containing a single derived {\tt type} definition comprises a class. 
The names of the module and the derived type correspond to each other (up to an extra underscore preceding the suffix in the module name). 
Derived types with data components have been available since Fortran 90; these data components are the {\em members} of the class.
But derived types with procedure components are new to Fortran 2003; these procedure components are the {\em methods} of the class. 
The language allows for two kinds of procedure components: procedure pointer variables, and type-bound procedures. 
We generally use the latter.
Type-bound procedures are declared in the {\tt contains} section of the derived type definition, while the subroutine or functions themselves are placed in the {\tt contains} section of the module.

Features new to Fortran 2003 also allow us to implement the principles of inheritance and polymorphism.
The new {\tt extends} keyword allows a new derived type (a child) to directly inherit the members and methods of a previously defined derived type (the parent), so that these components can be referenced directly as if they had been declared in the new child class itself. 
The new derived type can add additional members and methods to those already available through the parent type.
The child also can modify methods of the parent in two ways.
First, with the new {\tt generic} keyword it can {\em overload} a type-bound procedure name that is resolved to one of multiple procedures depending on the argument list. 
Second, it can {\em override} a type-bound procedure with a new subroutine or function exhibiting different behavior than the one in the parent class that it replaces.
We also note the new capability of defining a derived type as {\tt abstract}, and that includes {\tt deferred} methods for which only an interface---and not the subroutine or function itself---are provided.
Objects of an abstract type cannot be instantiated (i.e. variables of an abstract type cannot be declared); only objects of an extension of the abstract type, with type-bound procedures that override the deferred methods, can be instantiated. 

These mechanisms---which prove so useful in developing solvers and physics capabilities for GenASiS---are illustrated in simple cases in Listings~\ref{lst:Module_SWD_Outline} and \ref{lst:Module_WD_Outline} from the {\tt SineWaveDerivative} example, and Listings~\ref{lst:Module_GG_Outline} and \ref{lst:Module_FG_Outline} from the {\tt GaussianGradient} example.
In each of these problems the derivative or gradient of only a single function is taken.
But in each pair of listings, the parent class is an abstract template that, while containing some 80\% of the lines of code, does not depend upon the particular functional form at all.
The child class in each pair of listings serves only to specify a particular functional form.
It is evident that it would be simple and easy to make several other extensions for other functional forms.
Of course, in these trivial examples, less sophisticated techniques (such as passing the name of a function as a subroutine argument) could achieve the same purpose;
but in more complex and realistic situations, trying to achieve similar generality in Fortran with pre-Fortran 2003 language features fast becomes difficult and unwieldy.
Language constructs enabling object-oriented programming provide a much better means of achieving the same sort of code reusability.
The point is that the principles of inheritance and polymorphism---exemplified here in the mechanisms of {\tt type} extension and method overriding---make it much easier to allow lower-level code to access higher-level code.
This language functionality greatly facilitates, for instance, the separation of the {\tt Mathematics} and {\tt Physics} divisions of GenASiS. 
For example, solvers for generic classes of equations can be written in {\tt Mathematics}, and then invoked later by a range of different systems whose details are specified in {\tt Physics}.
This tremendously enhances the ease and transparency with which one can develop (for example) versatile and widely-applicable solvers.

An object-oriented approach need not ruin code performance.
This is because the floating-point operations on which the bulk of computational effort is (hopefully) spent can be isolated in low-level kernels expressed in terms of elementary loops over operations on array variables of intrinsic data types, which the compiler can handle well.
Consider for example the main computational task in the {\tt SineWaveDerivative} example: computing the derivative of the wave.
The wave and derivative data are stored in a member of an object of {\tt VariableGroupForm}, in such a way (see ``VariableGroups'' in Section~\ref{subsubsec:VariableManagement}) that reference to that data may at first look intimidating---not only to humans, but perhaps also to compilers!
However the actual operations on that data look straightforward to both of these audiences:
\begin{lstlisting}
  do i = WD
    dfdx ( i ) = (  f ( i + 1 ) - f ( i - 1 )  ) / ( 2.0 * dx )
  end do
\end{lstlisting} 
Here {\tt f} and {\tt dfdx} are simply rank-one {\tt real} dummy arguments of a kernel subroutine. 
The actual arguments are columns of the variable group's rank-two data array (and are therefore contiguous in memory), but reference to them here is plenty simple enough for the compiler to be able to apply any optimizations that may be useful.
Once again, this example problem is trivial, this time in the sense that there is so little data on which to work; 
but this same principle, namely of isolating computational work in kernel subroutines with basic loops over operations on simple dummy arguments of intrinsic data types, continues to be useful when the amount of data is much greater.
Our computational kernels throughout GenASiS---including in the {\tt Calculus} classes that compute differences and gradients for the {\tt GaussianGradient} example---take this approach.

Subsequent papers in this series will document additions to {\tt Mathematics} and {\tt Physics} that will make GenASiS\ suitable for multiphysics simulations of core-collapse supernovae and other astrophysical problems.  
Development of additional solvers and physics are already underway.  
Capabilities for nonrelativistic hydrodynamics will be reported first \citep{Cardall2012}.
We have implemented magnetohydrodynamics (MHD) capabilities in a `draft version' of GenASiS\ \citep{Endeve2012b}, which has been used to study standing accretion shock instability (SASI)-driven magnetic field amplification \citep{Endeve2010,Endeve2012}.  
In the current version of GenASiS\ we plan to implement MHD capabilities based on the HLLD solver \citep{Miyoshi2005,Mignone2009}.  
For Newtonian gravity, a multilevel Poisson solver will use a distributed FFT solver \citep{Budiardja2011} for the coarsest level, in conjunction with finite-difference solves on individual levels, in a multigrid approach.  
Work on general relativistic gravity---a major undertaking---is also underway  in GenASiS\ \citep{Tsatsin2011}.  
The greatest challenge of all is neutrino transport.  
Building on our group's past experience \citep{Liebendorfer2004,Bruenn2009}, theoretical developments \citep{Cardall2003,Cardall2005b}, and initial forays into transport in the earliest version of GenASiS\ \citep{Cardall2005,Cardall2009}, we plan to deploy a multigrid approach here as well, using a multigroup Variable Eddington Tensor formulation (Endeve et al. 2012, in preparation) with closures of increasing sophistication, ultimately culminating in a full `Boltzmann solver.'  


\acknowledgments

This research was supported by the Office of Advanced Scientific Computing Research and the Office of Nuclear Physics, U.S. Department of Energy.  
This research used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory provided through the INCITE program.  


\bibliographystyle{apalike}
