rocm-systems/examples/lulesh/lulesh-init.cc

#include <math.h>
#if USE_MPI
#    include <mpi.h>
#endif
#if _OPENMP
#    include <omp.h>
#endif
#include "lulesh.h"
#include <cstdlib>
#include <limits.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

static KOKKOS_INLINE_FUNCTION Real_t
                              CalcElemVolume(const Real_t x0, const Real_t x1, const Real_t x2, const Real_t x3,
                                             const Real_t x4, const Real_t x5, const Real_t x6, const Real_t x7,
                                             const Real_t y0, const Real_t y1, const Real_t y2, const Real_t y3,
                                             const Real_t y4, const Real_t y5, const Real_t y6, const Real_t y7,
                                             const Real_t z0, const Real_t z1, const Real_t z2, const Real_t z3,
                                             const Real_t z4, const Real_t z5, const Real_t z6, const Real_t z7)
{
    Real_t twelveth = Real_t(1.0) / Real_t(12.0);

    Real_t dx61 = x6 - x1;
    Real_t dy61 = y6 - y1;
    Real_t dz61 = z6 - z1;

    Real_t dx70 = x7 - x0;
    Real_t dy70 = y7 - y0;
    Real_t dz70 = z7 - z0;

    Real_t dx63 = x6 - x3;
    Real_t dy63 = y6 - y3;
    Real_t dz63 = z6 - z3;

    Real_t dx20 = x2 - x0;
    Real_t dy20 = y2 - y0;
    Real_t dz20 = z2 - z0;

    Real_t dx50 = x5 - x0;
    Real_t dy50 = y5 - y0;
    Real_t dz50 = z5 - z0;

    Real_t dx64 = x6 - x4;
    Real_t dy64 = y6 - y4;
    Real_t dz64 = z6 - z4;

    Real_t dx31 = x3 - x1;
    Real_t dy31 = y3 - y1;
    Real_t dz31 = z3 - z1;

    Real_t dx72 = x7 - x2;
    Real_t dy72 = y7 - y2;
    Real_t dz72 = z7 - z2;

    Real_t dx43 = x4 - x3;
    Real_t dy43 = y4 - y3;
    Real_t dz43 = z4 - z3;

    Real_t dx57 = x5 - x7;
    Real_t dy57 = y5 - y7;
    Real_t dz57 = z5 - z7;

    Real_t dx14 = x1 - x4;
    Real_t dy14 = y1 - y4;
    Real_t dz14 = z1 - z4;

    Real_t dx25 = x2 - x5;
    Real_t dy25 = y2 - y5;
    Real_t dz25 = z2 - z5;

#define TRIPLE_PRODUCT(x1, y1, z1, x2, y2, z2, x3, y3, z3)                               \
    ((x1) * ((y2) * (z3) - (z2) * (y3)) + (x2) * ((z1) * (y3) - (y1) * (z3)) +           \
     (x3) * ((y1) * (z2) - (z1) * (y2)))

    Real_t volume = TRIPLE_PRODUCT(dx31 + dx72, dx63, dx20, dy31 + dy72, dy63, dy20,
                                   dz31 + dz72, dz63, dz20) +
                    TRIPLE_PRODUCT(dx43 + dx57, dx64, dx70, dy43 + dy57, dy64, dy70,
                                   dz43 + dz57, dz64, dz70) +
                    TRIPLE_PRODUCT(dx14 + dx25, dx61, dx50, dy14 + dy25, dy61, dy50,
                                   dz14 + dz25, dz61, dz50);

#undef TRIPLE_PRODUCT

    volume *= twelveth;

    return volume;
}

/******************************************/

KOKKOS_INLINE_FUNCTION
Real_t
CalcElemVolume(const Real_t x[8], const Real_t y[8], const Real_t z[8])
{
    return CalcElemVolume(x[0], x[1], x[2], x[3], x[4], x[5], x[6], x[7], y[0], y[1],
                          y[2], y[3], y[4], y[5], y[6], y[7], z[0], z[1], z[2], z[3],
                          z[4], z[5], z[6], z[7]);
}

/////////////////////////////////////////////////////////////////////
Domain::Domain(Int_t numRanks, Index_t colLoc, Index_t rowLoc, Index_t planeLoc,
               Index_t nx, int tp, int nr, int balance, Int_t cost)
: m_e_cut(Real_t(1.0e-7))
, m_p_cut(Real_t(1.0e-7))
, m_q_cut(Real_t(1.0e-7))
, m_v_cut(Real_t(1.0e-10))
, m_u_cut(Real_t(1.0e-7))
, m_hgcoef(Real_t(3.0))
, m_ss4o3(Real_t(4.0) / Real_t(3.0))
, m_qstop(Real_t(1.0e+12))
, m_monoq_max_slope(Real_t(1.0))
, m_monoq_limiter_mult(Real_t(2.0))
, m_qlc_monoq(Real_t(0.5))
, m_qqc_monoq(Real_t(2.0) / Real_t(3.0))
, m_qqc(Real_t(2.0))
, m_eosvmax(Real_t(1.0e+9))
, m_eosvmin(Real_t(1.0e-9))
, m_pmin(Real_t(0.))
, m_emin(Real_t(-1.0e+15))
, m_dvovmax(Real_t(0.1))
, m_refdens(Real_t(1.0))
,
    //
    // set pointers to (potentially) "new'd" arrays to null to
    // simplify deallocation.
    //
    m_regNumList(0)
, m_nodeElemStart(0)
, m_nodeElemCornerList(0)
, m_regElemSize(0)
, m_regElemlist(0)
#if USE_MPI
, commDataSend(0)
, commDataRecv(0)
#endif
{
    Index_t edgeElems = nx;
    Index_t edgeNodes = edgeElems + 1;
    this->cost()      = cost;

    m_tp       = tp;
    m_numRanks = numRanks;

    ///////////////////////////////
    //   Initialize Sedov Mesh
    ///////////////////////////////

    // construct a uniform box for this processor

    m_colLoc   = colLoc;
    m_rowLoc   = rowLoc;
    m_planeLoc = planeLoc;

    m_sizeX   = edgeElems;
    m_sizeY   = edgeElems;
    m_sizeZ   = edgeElems;
    m_numElem = edgeElems * edgeElems * edgeElems;

    m_numNode = edgeNodes * edgeNodes * edgeNodes;

    m_regNumList = Allocate<Index_t>(numElem());  // material indexset

    // Elem-centered
    AllocateElemPersistent(numElem());

    // Node-centered
    AllocateNodePersistent(numNode());

    SetupCommBuffers(edgeNodes);

    // Basic Field Initialization
    for(Index_t i = 0; i < numElem(); ++i)
    {
        e(i)  = Real_t(0.0);
        p(i)  = Real_t(0.0);
        q(i)  = Real_t(0.0);
        ss(i) = Real_t(0.0);
    }

    // Note - v initializes to 1.0, not 0.0!
    for(Index_t i = 0; i < numElem(); ++i)
    {
        v(i) = Real_t(1.0);
    }

    for(Index_t i = 0; i < numNode(); ++i)
    {
        xd(i) = Real_t(0.0);
        yd(i) = Real_t(0.0);
        zd(i) = Real_t(0.0);
    }

    for(Index_t i = 0; i < numNode(); ++i)
    {
        xdd(i) = Real_t(0.0);
        ydd(i) = Real_t(0.0);
        zdd(i) = Real_t(0.0);
    }

    for(Index_t i = 0; i < numNode(); ++i)
    {
        nodalMass(i) = Real_t(0.0);
    }

    BuildMesh(nx, edgeNodes, edgeElems);

#if _OPENMP
    SetupThreadSupportStructures();
#else
    // These arrays are not used if we're not threaded
    m_nodeElemStart      = NULL;
    m_nodeElemCornerList = NULL;
#endif

    // Setup region index sets. For now, these are constant sized
    // throughout the run, but could be changed every cycle to
    // simulate effects of ALE on the lagrange solver
    CreateRegionIndexSets(nr, balance);

    // Setup symmetry nodesets
    SetupSymmetryPlanes(edgeNodes);

    // Setup element connectivities
    SetupElementConnectivities(edgeElems);

    // Setup symmetry planes and free surface boundary arrays
    SetupBoundaryConditions(edgeElems);

    // Setup defaults

    // These can be changed (requires recompile) if you want to run
    // with a fixed timestep, or to a different end time, but it's
    // probably easier/better to just run a fixed number of timesteps
    // using the -i flag in 2.x

    dtfixed()  = Real_t(-1.0e-6);  // Negative means use courant condition
    stoptime() = Real_t(1.0e-2);   // *Real_t(edgeElems*tp/45.0) ;

    // Initial conditions
    deltatimemultlb() = Real_t(1.1);
    deltatimemultub() = Real_t(1.2);
    dtcourant()       = Real_t(1.0e+20);
    dthydro()         = Real_t(1.0e+20);
    dtmax()           = Real_t(1.0e-2);
    time()            = Real_t(0.);
    cycle()           = Int_t(0);

    // initialize field data
    for(Index_t i = 0; i < numElem(); ++i)
    {
        Real_t   x_local[8], y_local[8], z_local[8];
        Index_t* elemToNode = nodelist(i);
        for(Index_t lnode = 0; lnode < 8; ++lnode)
        {
            Index_t gnode  = elemToNode[lnode];
            x_local[lnode] = x(gnode);
            y_local[lnode] = y(gnode);
            z_local[lnode] = z(gnode);
        }

        // volume calculations
        Real_t volume = CalcElemVolume(x_local, y_local, z_local);
        volo(i)       = volume;
        elemMass(i)   = volume;
        for(Index_t j = 0; j < 8; ++j)
        {
            Index_t idx = elemToNode[j];
            nodalMass(idx) += volume / Real_t(8.0);
        }
    }

    // deposit initial energy
    // An energy of 3.948746e+7 is correct for a problem with
    // 45 zones along a side - we need to scale it
    const Real_t ebase = Real_t(3.948746e+7);
    Real_t       scale = (nx * m_tp) / Real_t(45.0);
    Real_t       einit = ebase * scale * scale * scale;
    if(m_rowLoc + m_colLoc + m_planeLoc == 0)
    {
        // Dump into the first zone (which we know is in the corner)
        // of the domain that sits at the origin
        e(0) = einit;
    }
    // set initial deltatime base on analytic CFL calculation
    deltatime() = (Real_t(.5) * cbrt(volo(0))) / sqrt(Real_t(2.0) * einit);

}  // End constructor

////////////////////////////////////////////////////////////////////////////////
Domain::~Domain()
{
    /*   Release(&m_regNumList);
       Release(&m_nodeElemStart);
       Release(&m_nodeElemCornerList);
       Release(&m_regElemSize);
       for (Index_t i=0 ; i<numReg() ; ++i) {
         Release(&m_regElemlist[i]);
       }
       Release(&m_regElemlist);

    #if USE_MPI
       Release(&commDataSend);
       Release(&commDataRecv);
    #endif
    */
}  // End destructor

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

////////////////////////////////////////////////////////////////////////////////
void
Domain::BuildMesh(Int_t nx, Int_t edgeNodes, Int_t edgeElems)
{
    Index_t meshEdgeElems = m_tp * nx;

    // initialize nodal coordinates
    Index_t nidx = 0;
    Real_t  tz   = Real_t(1.125) * Real_t(m_planeLoc * nx) / Real_t(meshEdgeElems);
    for(Index_t plane = 0; plane < edgeNodes; ++plane)
    {
        Real_t ty = Real_t(1.125) * Real_t(m_rowLoc * nx) / Real_t(meshEdgeElems);
        for(Index_t row = 0; row < edgeNodes; ++row)
        {
            Real_t tx = Real_t(1.125) * Real_t(m_colLoc * nx) / Real_t(meshEdgeElems);
            for(Index_t col = 0; col < edgeNodes; ++col)
            {
                x(nidx) = tx;
                y(nidx) = ty;
                z(nidx) = tz;
                ++nidx;
                // tx += ds ; // may accumulate roundoff...
                tx = Real_t(1.125) * Real_t(m_colLoc * nx + col + 1) /
                     Real_t(meshEdgeElems);
            }
            // ty += ds ;  // may accumulate roundoff...
            ty = Real_t(1.125) * Real_t(m_rowLoc * nx + row + 1) / Real_t(meshEdgeElems);
        }
        // tz += ds ;  // may accumulate roundoff...
        tz = Real_t(1.125) * Real_t(m_planeLoc * nx + plane + 1) / Real_t(meshEdgeElems);
    }

    // embed hexehedral elements in nodal point lattice
    Index_t zidx = 0;
    nidx         = 0;
    for(Index_t plane = 0; plane < edgeElems; ++plane)
    {
        for(Index_t row = 0; row < edgeElems; ++row)
        {
            for(Index_t col = 0; col < edgeElems; ++col)
            {
                Index_t* localNode = nodelist(zidx);
                localNode[0]       = nidx;
                localNode[1]       = nidx + 1;
                localNode[2]       = nidx + edgeNodes + 1;
                localNode[3]       = nidx + edgeNodes;
                localNode[4]       = nidx + edgeNodes * edgeNodes;
                localNode[5]       = nidx + edgeNodes * edgeNodes + 1;
                localNode[6]       = nidx + edgeNodes * edgeNodes + edgeNodes + 1;
                localNode[7]       = nidx + edgeNodes * edgeNodes + edgeNodes;
                ++zidx;
                ++nidx;
            }
            ++nidx;
        }
        nidx += edgeNodes;
    }
}

////////////////////////////////////////////////////////////////////////////////
void
Domain::SetupThreadSupportStructures()
{
    // set up node-centered indexing of elements
    Index_t* nodeElemCount = Allocate<Index_t>(numNode());

    for(Index_t i = 0; i < numNode(); ++i)
    {
        nodeElemCount[i] = 0;
    }

    for(Index_t i = 0; i < numElem(); ++i)
    {
        Index_t* nl = nodelist(i);
        for(Index_t j = 0; j < 8; ++j)
        {
            ++(nodeElemCount[nl[j]]);
        }
    }

    m_nodeElemStart = Allocate<Index_t>(numNode() + 1);

    m_nodeElemStart[0] = 0;

    for(Index_t i = 1; i <= numNode(); ++i)
    {
        m_nodeElemStart[i] = m_nodeElemStart[i - 1] + nodeElemCount[i - 1];
    }

    m_nodeElemCornerList = Allocate<Index_t>(m_nodeElemStart[numNode()]);

    for(Index_t i = 0; i < numNode(); ++i)
    {
        nodeElemCount[i] = 0;
    }

    for(Index_t i = 0; i < numElem(); ++i)
    {
        Index_t* nl = nodelist(i);
        for(Index_t j = 0; j < 8; ++j)
        {
            Index_t m                    = nl[j];
            Index_t k                    = i * 8 + j;
            Index_t offset               = m_nodeElemStart[m] + nodeElemCount[m];
            m_nodeElemCornerList[offset] = k;
            ++(nodeElemCount[m]);
        }
    }

    Index_t clSize = m_nodeElemStart[numNode()];
    for(Index_t i = 0; i < clSize; ++i)
    {
        Index_t clv = m_nodeElemCornerList[i];
        if((clv < 0) || (clv > numElem() * 8))
        {
            fprintf(
                stderr,
                "AllocateNodeElemIndexes(): nodeElemCornerList entry out of range!\n");
#if USE_MPI
            MPI_Abort(MPI_COMM_WORLD, -1);
#else
            exit(-1);
#endif
        }
    }

    Release<Index_t>(&nodeElemCount);
}

////////////////////////////////////////////////////////////////////////////////
void
Domain::SetupCommBuffers(Int_t edgeNodes)
{
    // allocate a buffer large enough for nodal ghost data
    Index_t maxEdgeSize = MAX(this->sizeX(), MAX(this->sizeY(), this->sizeZ())) + 1;
    m_maxPlaneSize      = CACHE_ALIGN_REAL(maxEdgeSize * maxEdgeSize);
    m_maxEdgeSize       = CACHE_ALIGN_REAL(maxEdgeSize);

    // assume communication to 6 neighbors by default
    m_rowMin   = (m_rowLoc == 0) ? 0 : 1;
    m_rowMax   = (m_rowLoc == m_tp - 1) ? 0 : 1;
    m_colMin   = (m_colLoc == 0) ? 0 : 1;
    m_colMax   = (m_colLoc == m_tp - 1) ? 0 : 1;
    m_planeMin = (m_planeLoc == 0) ? 0 : 1;
    m_planeMax = (m_planeLoc == m_tp - 1) ? 0 : 1;

#if USE_MPI
    // account for face communication
    Index_t comBufSize =
        (m_rowMin + m_rowMax + m_colMin + m_colMax + m_planeMin + m_planeMax) *
        m_maxPlaneSize * MAX_FIELDS_PER_MPI_COMM;

    // account for edge communication
    comBufSize +=
        ((m_rowMin & m_colMin) + (m_rowMin & m_planeMin) + (m_colMin & m_planeMin) +
         (m_rowMax & m_colMax) + (m_rowMax & m_planeMax) + (m_colMax & m_planeMax) +
         (m_rowMax & m_colMin) + (m_rowMin & m_planeMax) + (m_colMin & m_planeMax) +
         (m_rowMin & m_colMax) + (m_rowMax & m_planeMin) + (m_colMax & m_planeMin)) *
        m_maxEdgeSize * MAX_FIELDS_PER_MPI_COMM;

    // account for corner communication
    // factor of 16 is so each buffer has its own cache line
    comBufSize +=
        ((m_rowMin & m_colMin & m_planeMin) + (m_rowMin & m_colMin & m_planeMax) +
         (m_rowMin & m_colMax & m_planeMin) + (m_rowMin & m_colMax & m_planeMax) +
         (m_rowMax & m_colMin & m_planeMin) + (m_rowMax & m_colMin & m_planeMax) +
         (m_rowMax & m_colMax & m_planeMin) + (m_rowMax & m_colMax & m_planeMax)) *
        CACHE_COHERENCE_PAD_REAL;

    this->commDataSend = Allocate<Real_t>(comBufSize);
    this->commDataRecv = Allocate<Real_t>(comBufSize);
    // prevent floating point exceptions
    memset(this->commDataSend, 0, comBufSize * sizeof(Real_t));
    memset(this->commDataRecv, 0, comBufSize * sizeof(Real_t));
#endif

    // Boundary nodesets
    if(m_colLoc == 0)
        m_symmX.resize(edgeNodes * edgeNodes);
    if(m_rowLoc == 0)
        m_symmY.resize(edgeNodes * edgeNodes);
    if(m_planeLoc == 0)
        m_symmZ.resize(edgeNodes * edgeNodes);
}

////////////////////////////////////////////////////////////////////////////////
void
Domain::CreateRegionIndexSets(Int_t nr, Int_t balance)
{
#if USE_MPI
    Index_t myRank;
    MPI_Comm_rank(MPI_COMM_WORLD, &myRank);
    srand(myRank);
#else
    srand(0);
    Index_t myRank = 0;
#endif
    this->numReg()    = nr;
    m_regElemSize     = Allocate<Index_t>(numReg());
    m_regElemlist     = Allocate<Index_t*>(numReg());
    Index_t nextIndex = 0;
    // if we only have one region just fill it
    // Fill out the regNumList with material numbers, which are always
    // the region index plus one
    if(numReg() == 1)
    {
        while(nextIndex < numElem())
        {
            this->regNumList(nextIndex) = 1;
            nextIndex++;
        }
        regElemSize(0) = 0;
    }
    // If we have more than one region distribute the elements.
    else
    {
        Int_t   regionNum;
        Int_t   regionVar;
        Int_t   lastReg = -1;
        Int_t   binSize;
        Index_t elements;
        Index_t runto           = 0;
        Int_t   costDenominator = 0;
        Int_t*  regBinEnd       = Allocate<Int_t>(numReg());
        // Determine the relative weights of all the regions.  This is based off the -b
        // flag.  Balance is the value passed into b.
        for(Index_t i = 0; i < numReg(); ++i)
        {
            regElemSize(i) = 0;
            costDenominator += pow((i + 1), balance);  // Total sum of all regions weights
            regBinEnd[i] =
                costDenominator;  // Chance of hitting a given region is (regBinEnd[i] -
                                  // regBinEdn[i-1])/costDenominator
        }
        // Until all elements are assigned
        while(nextIndex < numElem())
        {
            // pick the region
            regionVar = rand() % costDenominator;
            Index_t i = 0;
            while(regionVar >= regBinEnd[i])
                i++;
            // rotate the regions based on MPI rank.  Rotation is Rank % NumRegions this
            // makes each domain have a different region with the highest representation
            regionNum = ((i + myRank) % numReg()) + 1;
            // make sure we don't pick the same region twice in a row
            while(regionNum == lastReg)
            {
                regionVar = rand() % costDenominator;
                i         = 0;
                while(regionVar >= regBinEnd[i])
                    i++;
                regionNum = ((i + myRank) % numReg()) + 1;
            }
            // Pick the bin size of the region and determine the number of elements.
            binSize = rand() % 1000;
            if(binSize < 773)
            {
                elements = rand() % 15 + 1;
            }
            else if(binSize < 937)
            {
                elements = rand() % 16 + 16;
            }
            else if(binSize < 970)
            {
                elements = rand() % 32 + 32;
            }
            else if(binSize < 974)
            {
                elements = rand() % 64 + 64;
            }
            else if(binSize < 978)
            {
                elements = rand() % 128 + 128;
            }
            else if(binSize < 981)
            {
                elements = rand() % 256 + 256;
            }
            else
                elements = rand() % 1537 + 512;
            runto = elements + nextIndex;
            // Store the elements.  If we hit the end before we run out of elements then
            // just stop.
            while(nextIndex < runto && nextIndex < numElem())
            {
                this->regNumList(nextIndex) = regionNum;
                nextIndex++;
            }
            lastReg = regionNum;
        }
    }
    // Convert regNumList to region index sets
    // First, count size of each region
    for(Index_t i = 0; i < numElem(); ++i)
    {
        int r = this->regNumList(i) - 1;  // region index == regnum-1
        regElemSize(r)++;
    }
    // Second, allocate each region index set
    for(Index_t i = 0; i < numReg(); ++i)
    {
        m_regElemlist[i] = Allocate<Int_t>(regElemSize(i));
        regElemSize(i)   = 0;
    }
    // Third, fill index sets
    for(Index_t i = 0; i < numElem(); ++i)
    {
        Index_t r              = regNumList(i) - 1;  // region index == regnum-1
        Index_t regndx         = regElemSize(r)++;   // Note increment
        regElemlist(r, regndx) = i;
    }
}

/////////////////////////////////////////////////////////////
void
Domain::SetupSymmetryPlanes(Int_t edgeNodes)
{
    Index_t nidx = 0;
    for(Index_t i = 0; i < edgeNodes; ++i)
    {
        Index_t planeInc = i * edgeNodes * edgeNodes;
        Index_t rowInc   = i * edgeNodes;
        for(Index_t j = 0; j < edgeNodes; ++j)
        {
            if(m_planeLoc == 0)
            {
                m_symmZ[nidx] = rowInc + j;
            }
            if(m_rowLoc == 0)
            {
                m_symmY[nidx] = planeInc + j;
            }
            if(m_colLoc == 0)
            {
                m_symmX[nidx] = planeInc + j * edgeNodes;
            }
            ++nidx;
        }
    }
}

/////////////////////////////////////////////////////////////
void
Domain::SetupElementConnectivities(Int_t edgeElems)
{
    lxim(0) = 0;
    for(Index_t i = 1; i < numElem(); ++i)
    {
        lxim(i)     = i - 1;
        lxip(i - 1) = i;
    }
    lxip(numElem() - 1) = numElem() - 1;

    for(Index_t i = 0; i < edgeElems; ++i)
    {
        letam(i)                         = i;
        letap(numElem() - edgeElems + i) = numElem() - edgeElems + i;
    }
    for(Index_t i = edgeElems; i < numElem(); ++i)
    {
        letam(i)             = i - edgeElems;
        letap(i - edgeElems) = i;
    }

    for(Index_t i = 0; i < edgeElems * edgeElems; ++i)
    {
        lzetam(i) = i;
        lzetap(numElem() - edgeElems * edgeElems + i) =
            numElem() - edgeElems * edgeElems + i;
    }
    for(Index_t i = edgeElems * edgeElems; i < numElem(); ++i)
    {
        lzetam(i)                         = i - edgeElems * edgeElems;
        lzetap(i - edgeElems * edgeElems) = i;
    }
}

/////////////////////////////////////////////////////////////
void
Domain::SetupBoundaryConditions(Int_t edgeElems)
{
    Index_t ghostIdx[6];  // offsets to ghost locations

    // set up boundary condition information
    for(Index_t i = 0; i < numElem(); ++i)
    {
        elemBC(i) = Int_t(0);
    }

    for(Index_t i = 0; i < 6; ++i)
    {
        ghostIdx[i] = INT_MIN;
    }

    Int_t pidx = numElem();
    if(m_planeMin != 0)
    {
        ghostIdx[0] = pidx;
        pidx += sizeX() * sizeY();
    }

    if(m_planeMax != 0)
    {
        ghostIdx[1] = pidx;
        pidx += sizeX() * sizeY();
    }

    if(m_rowMin != 0)
    {
        ghostIdx[2] = pidx;
        pidx += sizeX() * sizeZ();
    }

    if(m_rowMax != 0)
    {
        ghostIdx[3] = pidx;
        pidx += sizeX() * sizeZ();
    }

    if(m_colMin != 0)
    {
        ghostIdx[4] = pidx;
        pidx += sizeY() * sizeZ();
    }

    if(m_colMax != 0)
    {
        ghostIdx[5] = pidx;
    }

    // symmetry plane or free surface BCs
    for(Index_t i = 0; i < edgeElems; ++i)
    {
        Index_t planeInc = i * edgeElems * edgeElems;
        Index_t rowInc   = i * edgeElems;
        for(Index_t j = 0; j < edgeElems; ++j)
        {
            if(m_planeLoc == 0)
            {
                elemBC(rowInc + j) |= ZETA_M_SYMM;
            }
            else
            {
                elemBC(rowInc + j) |= ZETA_M_COMM;
                lzetam(rowInc + j) = ghostIdx[0] + rowInc + j;
            }

            if(m_planeLoc == m_tp - 1)
            {
                elemBC(rowInc + j + numElem() - edgeElems * edgeElems) |= ZETA_P_FREE;
            }
            else
            {
                elemBC(rowInc + j + numElem() - edgeElems * edgeElems) |= ZETA_P_COMM;
                lzetap(rowInc + j + numElem() - edgeElems * edgeElems) =
                    ghostIdx[1] + rowInc + j;
            }

            if(m_rowLoc == 0)
            {
                elemBC(planeInc + j) |= ETA_M_SYMM;
            }
            else
            {
                elemBC(planeInc + j) |= ETA_M_COMM;
                letam(planeInc + j) = ghostIdx[2] + rowInc + j;
            }

            if(m_rowLoc == m_tp - 1)
            {
                elemBC(planeInc + j + edgeElems * edgeElems - edgeElems) |= ETA_P_FREE;
            }
            else
            {
                elemBC(planeInc + j + edgeElems * edgeElems - edgeElems) |= ETA_P_COMM;
                letap(planeInc + j + edgeElems * edgeElems - edgeElems) =
                    ghostIdx[3] + rowInc + j;
            }

            if(m_colLoc == 0)
            {
                elemBC(planeInc + j * edgeElems) |= XI_M_SYMM;
            }
            else
            {
                elemBC(planeInc + j * edgeElems) |= XI_M_COMM;
                lxim(planeInc + j * edgeElems) = ghostIdx[4] + rowInc + j;
            }

            if(m_colLoc == m_tp - 1)
            {
                elemBC(planeInc + j * edgeElems + edgeElems - 1) |= XI_P_FREE;
            }
            else
            {
                elemBC(planeInc + j * edgeElems + edgeElems - 1) |= XI_P_COMM;
                lxip(planeInc + j * edgeElems + edgeElems - 1) = ghostIdx[5] + rowInc + j;
            }
        }
    }
}

///////////////////////////////////////////////////////////////////////////
void
InitMeshDecomp(Int_t numRanks, Int_t myRank, Int_t* col, Int_t* row, Int_t* plane,
               Int_t* side)
{
    Int_t testProcs;
    Int_t dx, dy, dz;
    Int_t myDom;

    // Assume cube processor layout for now
    testProcs = Int_t(cbrt(Real_t(numRanks)) + 0.5);
    if(testProcs * testProcs * testProcs != numRanks)
    {
        printf("Num processors must be a cube of an integer (1, 8, 27, ...)\n");
#if USE_MPI
        MPI_Abort(MPI_COMM_WORLD, -1);
#else
        exit(-1);
#endif
    }
    if(sizeof(Real_t) != 4 && sizeof(Real_t) != 8)
    {
        printf("MPI operations only support float and double right now...\n");
#if USE_MPI
        MPI_Abort(MPI_COMM_WORLD, -1);
#else
        exit(-1);
#endif
    }
    if(MAX_FIELDS_PER_MPI_COMM > CACHE_COHERENCE_PAD_REAL)
    {
        printf("corner element comm buffers too small.  Fix code.\n");
#if USE_MPI
        MPI_Abort(MPI_COMM_WORLD, -1);
#else
        exit(-1);
#endif
    }

    dx = testProcs;
    dy = testProcs;
    dz = testProcs;

    // temporary test
    if(dx * dy * dz != numRanks)
    {
        printf("error -- must have as many domains as procs\n");
#if USE_MPI
        MPI_Abort(MPI_COMM_WORLD, -1);
#else
        exit(-1);
#endif
    }
    Int_t remainder = dx * dy * dz % numRanks;
    if(myRank < remainder)
    {
        myDom = myRank * (1 + (dx * dy * dz / numRanks));
    }
    else
    {
        myDom = remainder * (1 + (dx * dy * dz / numRanks)) +
                (myRank - remainder) * (dx * dy * dz / numRanks);
    }

    *col   = myDom % dx;
    *row   = (myDom / dx) % dy;
    *plane = myDom / (dx * dy);
    *side  = testProcs;

    return;
}