grid/gridDriver.h

Go to the documentation of this file.

/**
 *  @file      gridDriver.h
 *  @brief     Implements a spatially cell-centered uniform grid.
 *  @author    Mikica Kocic
 *  @copyright GNU General Public License (GPLv3).
 *
 *  @page fdm Finite Difference Approximation
 *
 *  At the core of finite difference approximation is a discretization of the
 *  spacetime, or a numerical grid (shown below).\ Then, a function `f(t,r)` is
 *  represented by values at a discrete set of points `(m,n)`.
 *
 *  <pre>
 *                p in M          <-->     (t,r)      <-->        (m,n)
 *         (point in a manifold)        (in a chart)        (discretized point)
 *                                           |                      |
 *                                           V                      V
 *                                         f(t,r)                 f(m,n)    </pre>
 *
 *  The points `(m,n)` are called the grid points or nodes. The grid points are chosen
 *  to reside at the center of grid cells (the other approach would be to associate the
 *  grid points with the vertices of the grid cells).&nbsp;
 *
 *  The data structure GridPoint contains the set of fields `{f,...}` that are kept
 *  in the memory for the point `(m,n)`. The numerical grid is an array of grid points
 *  `grid[m][n]` where a field `f` is directly accessible as `grid[m][n].f`, or using
 *  a macro as `f(m,n)`.&nbsp;
 *
 *  The grid consists of the *interior points*, and the *virtual points* at the boundary.
 *  The virtual points are used to impose the boundary conditions.&nbsp;
 *
 *  A cell-centered uniform grid with `nLen` grid points in each grid row is shown below
 *  (see also @cite Baumgarte:2010nr, p.192).
 *
 *  <pre>
 *         left virtual                grid cell        grid point     right virtual
 *            point(s)                 .---^---.            |             point(s)
 *      .        |               n-1  '    n    '  n+1      |                |
 *      :   +----|----+-------+-------+---------+-------+---|---+-------+----|----+
 *     m+1  |    O    |       |       | (m+1,n) |       |   X   |       |    O    |
 *          +---------+-------+-------.=========.-------+-------+-------+---------+
 *      m   |         |       |(m,n-1)|  (m,n)  |(m,n+1)|       |       |         |
 *          +---------+-------+-------'========='-------+-------+-------+---------+
 *     m-1  |         |       |       |         |       |       |       |         |
 *      .   +---------'-------+-------+---------+-------+-------+-------'---------+
 *      :        0    |   1       2      . . .                    nLen  |  nLen+1
 *                    ^                                                 ^
 *                   r_min                                            r_max       </pre>
 *
 *  The cell-centered uniform grid implies:
 *   @code
 *     r = r_min + ( n - 1/2 ) delta_r
 *     delta_r = ( r_max - r_min ) / nLen
 *   @endcode
 *   @li if `nGhost = 1`, `n` ranges between 1 and `nLen`;
 *       the boundary virtual points are `grid[m][0]` and `grid[m][nLen+1]`
 *
 *   @li `m` ranges between 1 and mLen, and wraps around;
 *       the slices `grid[mLen+1]`... are used for intermediate integration steps.
 *
 *  In parallel simulations, the grid is split across processes using MPI,
 *  allowing storage of large grid functions and parallel work on each part.
 *
 */

/////////////////////////////////////////////////////////////////////////////////////////
/** @defgroup g4 Grid driver                                                           */
/////////////////////////////////////////////////////////////////////////////////////////
                                                                                   /*@{*/
/** Uniform grid of grid functions.
 */
class UniformGrid
    : public MPIWorld
{
    /** Allocates memory for the numerical grid.
     *
     *  We allocate `mExtra` extra m grid rows for the intermediate steps.
     *  We also allocate 2*nGhost n-cells as the L/R spatial boundary.
     */
    void createGrid ()
    {
        size_t mDim = mLen + mExtra;
        size_t nDim = nLen + 2 * nGhost;

        grid = new GridPoint* [ mDim ];
        grid[0] = new GridPoint[ mDim * nDim ];

        for( size_t m = 0; m < mDim - 1; ++m ) {
            grid[m+1] = grid[m] + nDim;
        }

        slog << "    Allocated " << mDim << " * " << nDim
             << " * " << sizeof(GridPoint) / sizeof(Real)
             << " * Real" << 8 * sizeof(Real)
             << " = " << mDim * nDim * sizeof(GridPoint)
             << " bytes" << std::endl
             << "    Tracking " << GFCNT << " grid functions at a point"
             << std::endl << std::endl;
    }

    /** Releases memory allocated for the numerical grid.
     */
    void deleteGrid ()
    {
        if ( ! grid ) return;

        delete[] grid[0];
        delete[] grid;
        grid = NULL;
    }

protected:

    GridPoint** grid;      //!< The grid storage
                           //   Spatial discretization:
    Int         nLen;      //!< Number of cells in `r`-direction
    Int         nGhost;    //!< Number of virtual (ghost) cells on the `r`-boundary
    Real        delta_r;   //!< The grid spacing accross `r`
                           //   Temporal discretization:
    Int         mLen;      //!< Number of cached cells in `t`-direction
    Int         mExtra;    //!< Number of extra cells (for integration substeps)

public:

    GridPoint** get_grid () { return grid; }

    Int  get_nLen    () const { return nLen;    }
    Int  get_nGhost  () const { return nGhost;  }
    Int  get_mLen    () const { return mLen;    }
    Real get_delta_r () const { return delta_r; }

    Int  mpiSize     () const { return size (); }
    Int  mpiRank     () const { return rank (); }

    /** Constructs the uniform grid as specified in the parameter file.
     */
    UniformGrid( Parameters& params )
    {
        params.get( "grid.nLen",    nLen,     10  );
        params.get( "grid.nGhost",  nGhost,   10  );
        params.get( "grid.delta_r", delta_r,  0.1 );
        params.get( "grid.mLen",    mLen,     5   );
        params.get( "grid.mExtra",  mExtra,   9   );

        // Split the data equally among all the ranks in the MPI world
        nLen = nLen / size();

        slog << "Uniform Grid Driver:" << std::endl << std::endl
             << "    nLen = " << nLen << ",  nGhost = " << nGhost
             << ",  delta_r = " << delta_r
             << ",  mLen = " << mLen << ",  mExtra = " << mExtra << std::endl
             << "    Compile-time: CFDS_ORDER = " << CFDS_ORDER << std::endl;

        #if _OPENMP
             slog << "    Using OpenMP: #procs = "
                  << omp_get_num_procs() << ", max #threads = "
                  << omp_get_max_threads() << std::endl;
        #endif

        // Sanity check that nGhost > CFDS_ORDER/2
        //
        if( nGhost < CFDS_ORDER/2 ) {
            slog << "Error: Ghost region does not fit the FD order." << std::endl;
            exit( -1 );
        }

        createGrid ();
    }

    /** Exchange the boundaries at the given time slice.
     */
    void exchangeBoundaries( Int m )
    {
        MPIWorld::exchangeBoundaries(
            &grid[m][ 0                      ],  // Left ghost (received)
            &grid[m][ nGhost                 ],  // Left edge (sent)
            &grid[m][ nGhost + nLen - nGhost ],  // Right_edge (sent)
            &grid[m][ nGhost + nLen          ],  // Right ghost (received)
            nGhost * GFCNT                       // Exchanged size of Reals
        );
    }
};

/////////////////////////////////////////////////////////////////////////////////////////

/** Access the given grid function data (a utility macro for the grid users).
 */
#define GF(f,m,n)  grid[m][n][f]

/** `emitField` is a macro that emits a method encapsulating
 *  a grid function that access data in a grid at a point given by
 *  `m` (discretized t-coordinate) and `n` (discretized r-coordinate).
 *  For example, `emitField(r)` defines `r(m,n)`.
 */
#define emitField(gf) \
    inline Real& gf( Int m, Int n ) { return grid[m][n][fld::gf]; }

/** GridUser contains cached variables from the grid-driver.
 *  @note Finite differences are defined for the grid users.
 */
struct GridUser
{
    UniformGrid* gridDriver;  //!< Grid object
    GridPoint** grid;         //!< Actual grid data

    Int nLen;           //!< Number of cells in `r`-direction
    Int mLen;           //!< Number of cells in `t`-direction
    Int nGhost;         //!< Number of virtual (ghost) cells on the `r`-boundary

    Real delta_r;       //!< The grid spacing accross `r`
    Real inv_delta_r;   //!< `1 / delta_r`
    Real inv_delta_rr;  //!< `1 / delta_r^2`

    emitField(t)
    emitField(r)

    GridUser( UniformGrid& ug )
    {
        grid    = ug.get_grid ();
        nLen    = ug.get_nLen ();
        mLen    = ug.get_mLen ();
        nGhost  = ug.get_nGhost ();
        delta_r = ug.get_delta_r ();

        // Calculate various constants involving 1/delta_r
        //
        inv_delta_r  = 1 / ug.get_delta_r ();
        inv_delta_rr = 1 / ( ug.get_delta_r () * ug.get_delta_r () );
    }

    /** Cubic spline smoother of a grid function.
     */
    void cubicSplineShmooth( Int m, Int gf, Int lin2n, Int cub2n )
    {
        // The accuracy is very low at low r. Assume that a few first cells are linear.
        //
        if( lin2n > 0 )
        {
            Real dydx = delta_r * GF(gf,m,nGhost+lin2n) / r(m,nGhost+lin2n);
            for ( Int i = 0; i < lin2n; ++i ) {
                GF(gf,m,nGhost+i) = ( i + 0.5 ) * dydx;
            }
        }

        // A six-point cubic spline to smooth the region around r = 0.
        //
        if( cub2n > 0 )
        {
            static CubicSpline spline( 6 );

            Real pts[6][2] = {
                { r(m,nGhost           ),   GF(gf,m,nGhost           ) },
                { r(m,nGhost + 3       ),   GF(gf,m,nGhost + 3       ) },
                { r(m,nGhost + cub2n/2 ),   GF(gf,m,nGhost + cub2n/2 ) },
                { r(m,nGhost + cub2n-6 ),   GF(gf,m,nGhost + cub2n-6 ) },
                { r(m,nGhost + cub2n-3 ),   GF(gf,m,nGhost + cub2n-3 ) },
                { r(m,nGhost + cub2n   ),   GF(gf,m,nGhost + cub2n   ) }
            };
            spline.initialize( pts );

            for( Int i = 0; i < cub2n; ++i ) {
                GF(gf,m,nGhost+i) = spline( r(m,nGhost+i) );
            }
        }
    }
};
                                                                                   /*@}*/
/////////////////////////////////////////////////////////////////////////////////////////