*************** Software Design *************** Code structure ============== Castro is built upon the AMReX C++ framework. This provides high-level classes for managing an adaptive mesh refinement simulation, including the core data structures we will deal with. AMReX provides convenient data structures that allow for this workflow—high level objects in C++ that communicate with our computational kernels, with the data regions appearing as multidimensional arrays. Castro uses a structured-grid approach to hydrodynamics. We work with square/cubic zones that hold state variables (density, momentum, etc.) and compute the fluxes of these quantities through the interfaces of the zones (this is a finite-volume approach). Parallelization is achieved by domain decomposition. We divide our domain into many smaller boxes, and distributed these across processors. When information is needed at the boundaries of the boxes, messages are exchanged and this data is held in a perimeter of *ghost cells*. AMReX manages this decomposition and communication for us. Additionally, AMReX implements adaptive mesh refinement. In addition to the coarse decomposition of our domain into zones and boxes, we can refine rectangular regions by adding finer-gridded boxes on top of the coarser grid. We call the collection of boxes at the same resolution a *level*. Castro uses a hybrid MPI + OpenMP approach to parallelism. MPI is at used to communicate across nodes on a computer and OpenMP is used within a node, to loop over subregions of a box with different threads. The code structure in the Castro/ directory is: - ``Diagnostics/``: various analysis routines for specific problems - ``Docs/``: you’re reading this now! - ``Exec/``: various problem implementations, sorted by category: - ``gravity_tests/``: test problems that primarily exercise the gravity solver - ``hydro_tests/``: test problems of the hydrodynamics (with or without reactions) - ``radiation_tests/``: test problems that primarily exercise the radiation hydrodynamics solver - ``reacting_tests/``: test problems that primarily exercise the reactions (and hydro + reaction coupling) - ``scf_tests/``: problem setups that use the self-consistent field initialization - ``science/``: problem setups that were used for scientific investigations - ``unit_tests/``: test problems that exercise primarily a single module - ``external/``: if you are using git submodules, the Microphysics and AMReX git submodules will be in this directory. - ``paper/``: the JOSS paper source - ``Source/``: source code. In this main directory is all of the code. Sources are organized by topic as: - ``diffusion/`` : thermal diffusion code - ``driver/`` : the main driver, I/O, runtime parameter support - ``gravity/`` : self-gravity code - ``hydro/`` : the compressible hydrodynamics code - ``mhd/`` : the MHD solver code - ``particles/`` : support for particles - ``problems/`` : template code for implementing a problem - ``radiation/`` : the implicit radiation solve code - ``reactions/`` : nuclear reaction code - ``rotation/`` : rotating code - ``scf/`` : the self-consistent field initialization support - ``sdc/``: code specified for the true SDC method - ``sources/`` : hydrodynamics source terms support - ``Util/``: a catch-all for additional things you may need - ``ConvertCheckpoint/``: a tool to convert a checkpoint file to a larger domain - ... Major data structures ===================== The following data structures are the most commonly encountered when working in the C++ portions of Castro. This are all AMReX data-structures / classes. ``Amr`` ------- This is the main class that drives the whole simulation. This is the highest level in Castro. ``AmrLevel`` and Castro classes ------------------------------- An ``AmrLevel`` is a virtual base class provided by AMReX that stores all the state data on a single level in the AMR hierarchy and understands how to advance that data in time. The most important data managed by the ``AmrLevel`` is an array of ``StateData``, which holds the fluid quantities, etc., in the boxes that together make up the level. The ``Castro`` class is derived from the ``AmrLevel``. It provides the Castro-specific routines to evolve our system of equations. Like the ``AmrLevel``, there is one ``Castro`` object for each level in the AMR hierarchry. A lot of the member data in the ``Castro`` class are static member variables—this means that they are shared across all instances of the class. So, in this case, every level will have the same data. This is done, in particular, for the values of the runtime parameters, but also for the ``Gravity``, ``Diffusion``, and ``Radiation`` objects. This means that those objects cover all levels and are the same object in each instantiation of the ``Castro`` class. Floating point data ------------------- Floating point data in the C++ AMReX framework is declared as ``Real``. This is typedef to either ``float`` or ``double`` depending on the make variable ``PRECISION``. .. note :: single precision support in Castro is not yet complete. In particular, a lot of the supporting microphysics has not been updated. ``Box`` and ``FArrayBox`` ------------------------- A ``Box`` is simply a rectangular region in space. It does not hold data. In AMReX, an AMR level has a global index space, with :math:`(0,0,0)` being the lower left corner of the domain at that level, and :math:`(N_x-1, N_y-1, N_z-1)` being the upper right corner of the domain (for a domain of :math:`N_x \times N_y \times N_z` zones). The location of any ``Box`` at a level can be uniquely specified with respect to this global index space by giving the index of its lower-left and upper-right corners. :numref:`fig:soft:indexspace` shows an example of three boxes at the same level of refinement. AMReX provides other data structures that collect Boxes together, most importantly the ``BoxArray``. We generally do not use these directly, with the exception of the ``BoxArray`` ``grids``, which is defined as part of the ``AmrLevel`` class that ``Castro`` inherits. ``grids`` is used when building new ``MultiFabs`` to give the layout of the boxes at the current level. .. _fig:soft:indexspace: .. figure:: index_grid2.png :width: 4in Three boxes that comprise a single level. At this resolution, the domain is 20 :math:`\times` 18 zones. Note that the indexing in AMReX starts with :math:`0`. A ``FArrayBox`` or *FAB*, for *Fortran array box* is a data structure that contains a ``Box`` locating it in space, as well as a pointer to a data buffer. The real floating point data are stored as one-dimensional arrays in ``FArrayBox`` es. The associated ``Box`` can be used to reshape the 1D array into multi-dimensional arrays to be used by Fortran subroutines. The key part of the C++ AMReX data structures is that this data buffer can be cast a a DIM+1 dimensional array that we can easily fill in C++ kernels. .. note:: Castro is complied for a specific dimensionality. ``MultiFab`` ------------ At the highest abstraction level, we have the ``MultiFab`` (multiple FArrayBoxes). A ``MultiFab`` contains an array of ``Box`` es, including boxes owned by other processors for the purpose of communication, an array of MPI ranks specifying which MPI processor owns each ``Box``, and an array of pointers to ``FArrayBoxes`` owned by this MPI processor. .. note:: a ``MultiFab`` is a collection of the boxes that together make up a single level of data in the AMR hierarchy. A ``MultiFab`` can have multiple components (like density, temperature, ...) as well as a perimeter of ghost cells to exchange data with neighbors or implement boundary conditions (this is all reflected in the underlying ``FArrayBox``). Parallelization in AMReX is done by distributing the FABs across processors. Each processor knows which FABs are local to it. To loop over all the boxes local to a processor, an ``MFIter`` is used (more on this below). High-level operations exist on ``MultiFab`` s to add, subtract, multiply, etc., them together or with scalars, so you don’t need to write out loops over the data directly. In Castro, ``MultiFab`` s are one of the main data structures you will interact with in the C++ portions of the code. .. _soft:sec:statedata: ``StateData`` ------------- .. index:: StateData ``StateData`` is a class that essentially holds a pair of ``MultiFab`` s: one at the old time and one at the new time. AMReX knows how to interpolate in time between these states to get data at any intermediate point in time. The main data that we care about in Castro (the fluid state, gravitational potential, etc.) will be stored as ``StateData``. Essentially, data is made StateData in Castro if we need it to be stored in checkpoints / plotfiles, and/or we want it to be automatically interpolated when we refine. An ``AmrLevel`` stores an array of ``StateData`` (in a C++ array called ``state``). We index this array using integer keys (defined via an enum in ``Castro.H``). The state data is registered with AMReX in ``Castro_setup.cpp``. Note that each of the different ``StateData`` carried in the state array can have different numbers of components, ghost cells, boundary conditions, etc. This is the main reason we separate all this data into separate StateData objects collected together in an indexable array. The current ``StateData`` names Castro carries are: - ``State_Type`` : this is the ``NUM_STATE`` hydrodynamics components that make up the conserved hydrodynamics state (usually referred to as :math:`\Ub` in these notes. But note that this does not include the radiation energy density. We access this data using an AMReX ``Array4`` type which is of the form ``data(i,j,k,n)``, where ``n`` is the component. The integer keys used to index the components are defined in ``Source/driver/_variables`` (e.g., ``URHO``, ``UMX``, ``UMY``, ...) .. note:: regardless of dimensionality, we always carry around all three velocity components. The “out-of-plane” components will simply be advected, but we will allow rotation (in particular, the Coriolis force) to affect them. ``State_Type`` ``MultiFab`` s have no ghost cells. Note that the prediction of the hydrodynamic state to the interface will require 4 ghost cells. This accommodated by creating a separate MultiFab, ``Sborder`` that lives at the old-time level and has the necessary ghost cells. We will describe this more later. - ``Rad_Type`` : this stores the radiation energy density, commonly denoted :math:`E_r` in these notes. It has ``nGroups`` components—the number of energy groups used in the multigroup radiation hydrodynamics approximation. - ``PhiGrav_Type`` : this is simply the gravitational potential, usually denoted :math:`\Phi` in these notes. - ``Gravity_Type`` : this is the gravitational acceleration. There are always 3 components, regardless of the dimensionality (consistent with our choice of always carrying all 3 velocity components). - ``Source_Type`` : this holds the time-rate of change of the source terms, :math:`d\Sb/dt`, for each of the ``NUM_STATE`` ``State_Type`` variables. .. note:: we do not make use of the old-time quantity here. In fact, we never allocate the ``FArrayBox`` s for the old-time in the ``Source_Type`` ``StateData``, so there is not wasted memory. - ``Reactions_Type`` : this holds the data for the nuclear reactions. It has ``NumSpec+2`` components: the species creation rates (usually denoted :math:`\omegadot_k` in these notes), the specific energy generation rate (:math:`\dot{e}_\mathrm{nuc}`), and its density (:math:`\rho \dot{e}_\mathrm{nuc}`). These are stored as ``StateData`` so we have access to the reaction terms outside of advance, both for diagnostics (like flame speed estimation) and for reaction timestep limiting (this in particular needs the data stored in checkpoints for continuity of timestepping upon restart). - ``Mag_Type_x`` : this is defined for MHD and stores the face-centered (on x-faces) x-component of the magnetic field. - ``Mag_Type_y`` : this is defined for MHD and stores the face-centered (on y-faces) y-component of the magnetic field. - ``Mag_Type_z`` : this is defined for MHD and stores the face-centered (on z-faces) z-component of the magnetic field. - ``Simplified_SDC_React_Type`` : this is used with the SDC time-advancement algorithm. This stores the ``NQSRC`` terms that describe how the primitive variables change over the timestep due only to reactions. These are used when predicting the interface states of the primitive variables for the hydrodynamics portion of the algorithm. We access the ``MultiFab`` s that carry the data of interest by interacting with the ``StateData`` using one of these keys. For instance:: MultiFab& S_new = get_new_data(State_Type); gets a pointer to the ``MultiFab`` containing the hydrodynamics state data at the new time. ``MFIter`` ========== The process of looping over boxes at a given level of refinement and operating on their data is linked to how Castro achieves thread-level parallelism. The OpenMP approach in Castro has evolved considerably since the original paper was written, with the modern approach, called *tiling*, gearing up to meet the demands of many-core processors in the next-generation of supercomputers. Full details of iterating over boxes and calling compute kernels is given in the AMReX documentation here: https://amrex-codes.github.io/amrex/docs_html/Basics.html#mfiter-and-tiling Practical Details in Working with Tiling ---------------------------------------- With tiling, the OpenMP is now all at the loop over boxes and not in the computational kernels themselves. It is the responsibility of the coder to make sure that the routines within a tiled region are safe to use with OpenMP. In particular, note that: - tile boxes are non-overlapping - the union of tile boxes completely cover the valid region of the fab - Consider working with a node-centered MultiFab, ``ugdnv``, and a cell-centered ``MultiFab`` ``s``: - with ``mfi(s)``, the tiles are based on the cell-centered index space. If you have an :math:`8\times 8` box, then and 4 tiles, then your tiling boxes will range from :math:`0\rightarrow 3`, :math:`4\rightarrow 7`. - with ``mfi(ugdnv)``, the tiles are based on nodal indices, so your tiling boxes will range from :math:`0\rightarrow 3`, :math:`4\rightarrow 8`. - When updating routines to work with tiling, we need to understand the distinction between the index-space of the entire box (which corresponds to the memory layout) and the index-space of the tile. Boundaries: ``FillPatch`` and ``FillPatchIterator`` =================================================== AMReX calls the act of filling ghost cells a *fillpatch* operation. Boundaries between grids are of two types. The first we call “fine-fine”, which is two grids at the same level. The second type is "coarse-fine", which needs interpolation from the coarse grid to fill the fine grid ghost cells. Both of these are part of the fillpatch operation. Fine-fine fills are just a straight copy from “valid regions” to ghost cells. Coarse-fine fills are enabled because the ``StateData`` is not just arrays, they’re “State Data”, which means that the data knows how to interpolate itself (in an anthropomorphical sense). The type of interpolation to use is defined in ``Castro_setup.cpp``—search for ``cell_cons_interp``, for example—that’s “cell conservative interpolation”, i.e., the data is cell-based (as opposed to node-based or edge-based) and the interpolation is such that the average of the fine values created is equal to the coarse value from which they came. (This wouldn’t be the case with straight linear interpolation, for example.) Additionally, since ``StateData`` has an old and new timelevel, the fill patch operation can interpolate to an intermediate time. Examples -------- To illustrate the various ways we fill ghost cells and use the data, let’s consider the following scenarios: - *You have state data that was defined with no ghost cells. You want to create a new* ``MultiFab`` *containing a copy of that data with* ``NGROW`` *ghost cells.* This is the case with ``Sborder`` —the ``MultiFab`` of the hydrodynamic state that we use to kick-off the hydrodynamics advance. ``Sborder`` is declared in ``Castro.H`` simply as: .. code:: c++ Multifab Sborder; It is then allocated in ``Castro::initialize_do_advance()`` .. code:: c++ Sborder.define(grids, NUM_STATE, NUM_GROW, Fab_allocate); const Real prev_time = state[State_Type].prevTime(); expand_state(Sborder, prev_time, NUM_GROW); Note in the call to ``.define()``, we tell AMReX to already allocate the data regions for the ``FArrayBox`` s that are part of ``Sborder``. The actually filling of the ghost cells is done by ``Castro::expand_state()``: .. code:: c++ AmrLevel::FillPatch(*this, Sborder, NUM_GROW, prev_time, State_Type, 0, NUM_STATE); Here, we are filling the ng ghost cells of ``MultiFab`` ``Sborder`` at time prev_time. We are using the ``StateData`` that is part of the current ``Castro`` object that we are part of. Note: ``FillPatch`` takes an object reference as its first argument, which is the object that contains the relevant ``StateData`` —that is what the this pointer indicates. Finally, we are copying the ``State_Type`` data components 0 to ``NUM_STATE`` [1]_. The result of this operation is that ``Sborder`` will now have ``NUM_GROW`` ghost cells consistent with the ``State_Type`` data at the old time-level. - *You have state data that was defined with* ``NGROW`` *ghost cells. You want to ensure that the ghost cells are filled (including any physical boundaries) with valid data.* This is very similar to the procedure shown above. The main difference is that for the ``MultiFab`` that will be the target of the ghost cell filling, we pass in a reference to the ``StateData`` itself. The main thing you need to be careful of here, is that you need to ensure that the the time you are at is consistent with the ``StateData`` ’s time. Here’s an example from the radiation portion of the code ``MGFLDRadSolver.cpp``: .. code:: c++ Real time = castro->get_state_data(Rad_Type).curTime(); MultiFab& S_new = castro->get_new_data(State_Type); AmrLevel::FillPatch(*castro, S_new, ngrow, time, State_Type, 0, S_new.nComp(), 0); In this example, ``S_new`` is a pointer to the new-time-level ``State_Type`` ``MultiFab``. So this operation will use the ``State_Type`` data to fill its own ghost cells. we fill the ``ngrow`` ghost cells of the new-time-level ``State_Type`` data, for all the components. Note that in this example, because the ``StateData`` lives in the ``castro`` object and we are working from the ``Radiation`` object, we need to make reference to the current ``castro`` object pointer. If this were all done within the ``castro`` object, then the pointer will simply be ``this``, as we saw above. - *You have a* ``MultiFab`` *with some derived quantity. You want to fill its ghost cells.* ``MultiFabs`` have a ``FillBoundary()`` method that will fill all the ghost cells between boxes at the same level. It will not fill ghost cells at coarse-fine boundaries or at physical boundaries. - *You want to loop over the FABs in state data, filling ghost cells along the way* This is the job of the ``FillPatchIterator``—this iterator is used to loop over the grids and fill ghostcells. A key thing to keep in mind about the ``FillPatchIterator`` is that you operate on a copy of the data—the data is disconnected from the original source. If you want to update the data in the source, you need to explicitly copy it back. Also note: ``FillPatchIterator`` takes a ``MultiFab``, but this is not filled—this is only used to get the grid layout. Finally, the way the ``FillPatchIterator`` is implemented is that all the communication is done first, and then the iterating over boxes commences. For example, the loop that calls ``CA_UMDRV`` (all the hydrodynamics integration stuff) starts with:: for (FillPatchIterator fpi(*this, S_new, NUM_GROW, time, State_Type, strtComp, NUM_STATE); fpi.isValid(); ++fpi) { FArrayBox &state = fpi(); Box bx(fpi.validbox()); // work on the state FAB. The interior (valid) cells will // live between bx.loVect() and bx.hiVect() } Here the ``FillPatchIterator`` is the thing that distributes the grids over processors and makes parallel “just work”. This fills the single patch ``fpi`` , which has ``NUM_GROW`` ghost cells, with data of type ``State_Type`` at time ``time``, starting with component strtComp and including a total of ``NUM_STATE`` components. In general, one should never assume that ghostcells are valid, and instead do a fill patch operation when in doubt. Sometimes we will use a ``FillPatchIterator`` to fill the ghost cells into a ``MultiFab`` without an explicit look. This is done as:: FillPatchIterator fpi(*this,S_old,1,time,State_Type,0,NUM_STATE); MultiFab& state_old = fpi.get_mf(); In this operation, state_old points to the internal ``MultiFab`` in the ``FillPatchIterator``, by getting a reference to it as ``fpi.get_mf()``. This avoids a local copy. Note that in the examples above, we see that only ``StateData`` can fill physical boundaries (because these register how to fill the boundaries when they are defined). There are some advanced operations in AMReX that can get around this, but we do not use them in Castro. .. _soft:phys_bcs: Physical Boundaries ------------------- .. index:: boundary conditions Physical boundary conditions are specified by an integer index [2]_ in the ``inputs`` file, using the ``castro.lo_bc`` and ``castro.hi_bc`` runtime parameters. The generally supported boundary conditions are, their corresponding integer key, and the action they take for the normal velocity, transverse velocity, and generic scalar are shown in :numref:`table:castro:bcs`. The definition of the specific actions are: - ``INT_DIR``: data taken from other grids or interpolated - ``EXT_DIR``: data specified on EDGE (FACE) of bndry - ``HOEXTRAP``: higher order extrapolation to EDGE of bndry - ``FOEXTRAP``: first order extrapolation from last cell in interior - ``REFLECT_EVEN``: :math:`F(-n) = F(n)` true reflection from interior cells - ``REFLECT_ODD``: :math:`F(-n) = -F(n)` true reflection from interior cells The actual registration of a boundary condition action to a particular variable is done in ``Castro_setup.cpp``. At the top we define arrays such as ``scalar_bc``, ``norm_vel_bc``, etc, which say which kind of bc to use on which kind of physical boundary. Boundary conditions are set in functions like ``set_scalar_bc``, which uses the ``scalar_bc`` pre-defined arrays. We also specify the name of the routine that is responsible for filling the data there (e.g., ``hypfill``). These routines are discussed more below. If you want to specify a value at a function (like at an inflow boundary), then you choose an *inflow* boundary at that face of the domain. You then need to write the implementation code for this. There is a centralized hydrostatic boundary condition that is implemented this way—see :ref:`create:bcs`. .. _table:castro:bcs: .. table:: Physical boundary conditions supported in Castro. +-------------+-------------+-------------+--------------+--------------+ | **name** | **integer** | **normal | **transverse | **scalars** | | | | velocity** | velocity** | | +=============+=============+=============+==============+==============+ | interior | 0 | INT_DIR | INT_DIR | INT_DIR | +-------------+-------------+-------------+--------------+--------------+ | inflow | 1 | EXT_DIR | EXT_DIR | EXT_DIR | +-------------+-------------+-------------+--------------+--------------+ | outflow | 2 | FOEXTRAP | FOEXTRAP | FOEXTRAP | +-------------+-------------+-------------+--------------+--------------+ | symmetry | 3 | REFLECT_ODD | REFLECT_EVEN | REFLECT_EVEN | +-------------+-------------+-------------+--------------+--------------+ | slipwall | 4 | REFLECT_ODD | REFLECT_EVEN | REFLECT_EVEN | +-------------+-------------+-------------+--------------+--------------+ | noslipwall | 5 | REFLECT_ODD | REFLECT_EVEN | REFLECT_EVEN | +-------------+-------------+-------------+--------------+--------------+ ``FluxRegister`` ---------------- A ``FluxRegister`` holds face-centered data at the boundaries of a box. It is composed of a set of ``MultiFab`` s (one for each face, so 6 for 3D). A ``FluxRegister`` stores fluxes at coarse-fine interfaces, and isused for the flux-correction step. Other AMReX Concepts ==================== There are a large number of classes that help define the structure of the grids, metadata associate with the variables, etc. A good way to get a sense of these is to look at the ``.H`` files in the ``amrex/Src/`` directory. ``Geometry`` class ------------------ There is a ``Geometry`` object, ``geom`` for each level as part of the ``Castro`` object (this is inherited through ``AmrLevel``). ``ParmParse`` class ------------------- Error Estimators ---------------- ``Gravity`` class ================= There is a single ``Gravity`` object, ``gravity``, that is a static class member of the ``Castro`` object. This means that all levels refer to the same ``Gravity`` object. Within the ``Gravity`` object, there are pointers to the ``Amr`` object (as ``parent``), and all of the ``AmrLevels`` (as a ``PArray``, ``LevelData``). The ``gravity`` object gets the geometry information at each level through the parent ``Amr`` class. The main job of the ``gravity`` object is to provide the potential and gravitation acceleration for use in the hydrodynamic sources. Depending on the approximation used for gravity, this could mean calling the AMReX multigrid solvers to solve the Poisson equation. .. [1] for clarity and continuity in this documentation, some of the variable names have been changed compared to the actual code .. [2] the integer values are defined in ``BC_TYPES.H``