DgSlices< nDim, SysEqn > Class Template Reference

Determine all coordinates and alloc buffer size create a valid 2D grid which represents a slice from a 3D grid. More...

#include <dgcartesianslices.h>

Collaboration diagram for DgSlices< nDim, SysEqn >:

Public Member Functions
	DgSlices (DgCartesianSolver< nDim, SysEqn > &solver)
void	setProperties ()
	Read properties for all slices and create corresponding slice objects. More...
void	init ()
	Initialize all slice objects. More...
void	save (const MBool isFinalTimeStep)
	Call save method for all slice objects. More...
MBool	enabled () const

Private Member Functions
MPI_Comm	mpiComm () const
MInt	solverId () const
MInt	domainId () const
void	updateGridFile (DgSliceSeries &sliceSeries, const MString &prefix="")
void	init (DgSliceSeries &sliceSeries)
	Create 2D slice grid file and get all elements which belong to specific grid. More...
void	save (DgSliceSeries &sliceSeries, const MBool isFinalTimeStep)
	Collect data for all slice nodes and create solution files. Data is sorted by their slice hilbertId and therefore also written in chunks of hilbertIds. More...

Private Attributes
DgCartesianSolver< nDim, SysEqn > &	m_solver
MBool	m_enabled = false
std::vector< DgSliceSeries >	m_sliceSeries {}
std::vector< MInt >	m_sliceVarIds {}
	List of slice variables. More...
std::vector< MInt >	m_noSliceVars {}
	Number of variables for each slice variable. More...
std::vector< std::vector< MString > >	m_sliceVarNames {}
	List of variable names for each slice variable. More...

Static Private Attributes
static constexpr const MInt	m_dim = 3

Detailed Description

template<MInt nDim, class SysEqn>
class DgSlices< nDim, SysEqn >

Author

Marcus Wiens (marcus) marcu.nosp@m.s.wi.nosp@m.ens@r.nosp@m.wth-.nosp@m.aache.nosp@m.n.de

Date

06.01.2016

Definition at line 27 of file dgcartesianslices.h.

Constructor & Destructor Documentation

DgSlices()

template<MInt nDim, class SysEqn >

DgSlices< nDim, SysEqn >::DgSlices ( DgCartesianSolver< nDim, SysEqn > &  solver )

inlineexplicit

Definition at line 29 of file dgcartesianslices.h.

: m_solver(solver) {}

Member Function Documentation

domainId()

template<MInt nDim, class SysEqn >

MInt DgSlices< nDim, SysEqn >::domainId ( ) const

inlineprivate

Definition at line 38 of file dgcartesianslices.h.

{ return m_solver.domainId(); }

enabled()

template<MInt nDim, class SysEqn >

MBool DgSlices< nDim, SysEqn >::enabled ( ) const

inline

Definition at line 33 of file dgcartesianslices.h.

{ return m_enabled; }

init() [1/2]

template<MInt nDim, class SysEqn >

void DgSlices< nDim, SysEqn >::init

Author

Marcus Wiens (marcus) marcu.nosp@m.s.wi.nosp@m.ens@r.nosp@m.wth-.nosp@m.aache.nosp@m.n.de

Date

20.06.2016

Definition at line 339 of file dgcartesianslices.h.

                                  {
  TRACE();
  // Return if dg slice is deactivated
  if(!m_enabled) {
    return;
  }
 
  // Call init method for all slice objects
  for(MUint i = 0; i < m_sliceSeries.size(); i++) {
    init(m_sliceSeries[i]);
  }
}

init() [2/2]

template<MInt nDim, class SysEqn >

void DgSlices< nDim, SysEqn >::init ( DgSliceSeries &  sliceSeries )

private

   First the DG Elements, Number of Nodes, PolyDegree and Coorinates, which belong the
   current slice are determined from the slice cells. This is done in data chunks with the
   same hilbertId.

   After that a MPI Communicator is created for all domains of the current slice.

   In the last step, the HilbertId Information is used, to determine global offsets for file
   creation. Since the data in the grids and solution files is sorted by their hilbertId, it
   is necessary to determine this offset information. It is important to understand, that
   the number of hilbertIds is different on every domain and that they can be distributed in
   in any order, since a slice of a 3D hilbert curve is performed. To determine the offset
   correctly the local information of hilbertIds and number of nodes is exchanged  globally
   (globally in the context of slices), that every domain knows how many hilbertIds and
   nodes (with same hilbertId) are on each domain and from this the global offsets are
   calculated for the local data.

Author

Marcus Wiens (marcus) marcu.nosp@m.s.wi.nosp@m.ens@r.nosp@m.wth-.nosp@m.aache.nosp@m.n.de

Date

06.01.2016

Definition at line 376 of file dgcartesianslices.h.

                                                            {
  TRACE();
 
  // Reset variables, needed when DgSlices::init() is called multiple times, e.g., DLB enabled
  // Reset noLocalNodes
  sliceSeries.m_localNoNodes = 0;
 
  // Step 1: Determine DG Elements and number of Nodes
 
  // Create grid file and get slice cells from grid
  MIntScratchSpace cellIds(m_solver.grid().noInternalCells(), AT_, "cellIds");
  MIntScratchSpace hilbertInfo(m_solver.grid().noInternalCells() * 3, AT_, "hilbertIds");
  // Get cellIds for current slice and hilbertId informaiton. The cellIds are used to determine the
  // DG elements for the slice and the hilbertId information is needed for file creation.
  const MString gridFilePath = m_solver.outputDir() + sliceSeries.m_gridFileName;
  m_solver.grid().raw().createGridSlice(sliceSeries.m_axis, sliceSeries.m_intercept, gridFilePath, solverId(),
                                        &sliceSeries.m_noCellInSlice, &cellIds[0], &sliceSeries.m_noHilbertIds,
                                        &hilbertInfo[0]);
  sliceSeries.m_elementNodesPerHilbertIds.resize(sliceSeries.m_noHilbertIds);
  sliceSeries.m_elementIdsPerHilbertId.resize(sliceSeries.m_noHilbertIds);
 
  m_log << "Initializing DG Slice " << sliceSeries.m_fileNo << " ...";
  // Search for slice elements if slice cells are found on current domain
  if(sliceSeries.m_noCellInSlice > 0) {
    // Convert grid cell-ids into solver cell-ids (only required if there are multiple solvers)
    if(m_solver.grid().raw().treeb().noSolvers() > 1) {
      for(MInt i = 0; i < sliceSeries.m_noCellInSlice; i++) {
        cellIds[i] = m_solver.grid().tree().grid2solver(cellIds[i]);
      }
    }
 
    // Initalize polyDegree vector
    sliceSeries.m_polyDegs.resize(sliceSeries.m_noCellInSlice);
    std::fill_n(&sliceSeries.m_polyDegs[0], sliceSeries.m_noCellInSlice, m_solver.m_initPolyDeg);
 
    // Assume size of element variables
    const MInt noMaxElements = m_solver.m_elements.size();
    // Collectors for element data
    MIntScratchSpace elementIds(noMaxElements, AT_, "elementIds");
    MIntScratchSpace elementNodes(noMaxElements, AT_, "elementNodes");
    MIntScratchSpace elementOffset(noMaxElements, AT_, "elementOffset");
 
    // Total number of coordinates for all elements
    MInt noCoords = 0;
    // Assume max number of nodes
    const MInt noMaxNodes = pow(
        *std::max_element(&m_solver.m_elements.polyDeg(0), &m_solver.m_elements.polyDeg(0) + m_solver.m_elements.size())
            + 1,
        2);
    MFloatScratchSpace coordinates(noMaxElements * noMaxNodes * m_dim, AT_, "coordinates");
 
    // Holds coordinates for one element
    std::array<MFloat, m_dim> nodeCoord;
    // This loop searches for elements and node coordinates of the slice. The cells are checked in
    // chunks the same sliceHilbertId to determine the element node offset at the end of this part.
    MInt noElementsSlice = 0;
    for(MInt h = 0, cellSum = 0, noNodes = 0, prevNoNodes = 0, sumElements = 0; h < sliceSeries.m_noHilbertIds; h++) {
      for(MInt c = cellSum, e = 0; (c - cellSum) < hilbertInfo[h * 3 + 1]; e++) {
        // Check to identify cells which are not an element; These cells are
        // counted as dummy nodes in the solution file and have the value 0
 
        // skip slice which do not have a corresponding element
        if(cellIds[c] < m_solver.m_elements.cellId(e)) {
          // count dummy nodes to cells which are not elements
          noNodes += ipow(m_solver.m_initPolyDeg + 1, m_dim - 1);
          // Go to next cell after counting nodes
          c++;
          // Hold element until all dummy elements are counted
          e--;
          continue;
        }
        // save node coordinates if element is in slice
        if(cellIds[c] == m_solver.m_elements.cellId(e)) {
          // Determine number of nodes
          const MInt noNodes1D = m_solver.m_elements.polyDeg(e) + 1;
 
          // Tensor to hold node coordinates
          const MFloatTensor nodeCoordinates(&m_solver.m_elements.nodeCoordinates(e), noNodes1D, noNodes1D, noNodes1D,
                                             m_dim);
          // Loop over all nodes to save the coordinates
          for(MInt i = 0; i < noNodes1D; i++) {
            for(MInt j = 0; j < noNodes1D; j++) {
              if(sliceSeries.m_axis == "x") {
                nodeCoord[0] = sliceSeries.m_intercept;
                nodeCoord[1] = nodeCoordinates(0, i, j, 1);
                nodeCoord[2] = nodeCoordinates(0, i, j, 2);
              } else if(sliceSeries.m_axis == "y") {
                nodeCoord[0] = nodeCoordinates(i, 0, j, 0);
                nodeCoord[1] = sliceSeries.m_intercept;
                nodeCoord[2] = nodeCoordinates(i, 0, j, 2);
              } else if(sliceSeries.m_axis == "z") {
                nodeCoord[0] = nodeCoordinates(i, j, 0, 0);
                nodeCoord[1] = nodeCoordinates(i, j, 0, 1);
                nodeCoord[2] = sliceSeries.m_intercept;
              }
              // Save coordinate tuple
              std::copy_n(&nodeCoord[0], m_dim, &coordinates[noCoords]);
              noCoords += m_dim;
            }
          }
          elementIds[noElementsSlice] = e;
 
          // Otherwise use current offset for the current element
          elementOffset[noElementsSlice] = noNodes - prevNoNodes; // sliceSeries.m_localNoNodes;
 
          // Increase offset based on element polynomial degree
          const MInt noNodesXD = ipow(m_solver.m_elements.polyDeg(e) + 1, m_dim - 1);
          noNodes += noNodesXD;
 
          // Set number of nodes to use later and proceed with next element
          elementNodes[noElementsSlice] = noNodesXD;
          noElementsSlice++;
 
          // save polyDegree (it is easier to collect this data here)
          sliceSeries.m_polyDegs[c] = m_solver.m_elements.polyDeg(e);
 
          c++;
        }
      }
      // Save number of elements which have the same hilbertId
      sliceSeries.m_elementIdsPerHilbertId[h] = noElementsSlice - sumElements;
 
      // Save number of nodes per hilbert id and total number
      sliceSeries.m_elementNodesPerHilbertIds[h] = noNodes - prevNoNodes;
      sliceSeries.m_localNoNodes += sliceSeries.m_elementNodesPerHilbertIds[h];
      // Save number of evaluated cells, elements and nodes
      cellSum += hilbertInfo[h * 3 + 1];
      prevNoNodes = noNodes;
      sumElements = noElementsSlice;
    }
    // Resize variables
    sliceSeries.m_elementIds.resize(noElementsSlice);
    sliceSeries.m_elementNodes.resize(noElementsSlice);
    sliceSeries.m_elementOffset.resize(noElementsSlice);
    sliceSeries.m_coordinates.resize(noCoords);
    // Get values from scratches and save them
    std::copy_n(&elementIds[0], noElementsSlice, &sliceSeries.m_elementIds[0]);
    std::copy_n(&elementNodes[0], noElementsSlice, &sliceSeries.m_elementNodes[0]);
    std::copy_n(&elementOffset[0], noElementsSlice, &sliceSeries.m_elementOffset[0]);
    std::copy_n(&coordinates[0], noCoords, &sliceSeries.m_coordinates[0]);
  }
 
  // If only non element cells are found, make m_coordinates non zero length
  // to create a dummy states scratch
  if(sliceSeries.m_coordinates.empty()) {
    sliceSeries.m_coordinates.resize(1);
  }
 
  // Step 2: Create MPI Communicator
 
  using namespace maia;
  // Set MPI Comm for ranks that are part of slice
  // Create a new MPI Communicator
  ScratchSpace<MInt> sliceScratch(globalNoDomains(), AT_, "sliceScratch");
 
  // Get Rank if domain holds elements in slice
  MInt rank = -1;
  if(sliceSeries.m_noCellInSlice > 0) {
    MPI_Comm_rank(mpiComm(), &rank);
  }
 
  // Combine all ranks to get relevant ranks
  MPI_Allgather(&rank, 1, type_traits<MInt>::mpiType(), &sliceScratch[0], 1, type_traits<MInt>::mpiType(), mpiComm(),
                AT_, "rank", "sliceScratch[0]");
 
  // Check for relevant ranks and save them to create the new communicator
  const MInt noSliceDomains =
      std::count_if(sliceScratch.begin(), sliceScratch.end(), [](const MInt a) { return a != -1; });
  MIntScratchSpace sliceDomains(noSliceDomains, AT_, "sliceDomains");
  MInt position = 0;
  for(MInt i : sliceScratch) {
    if(i != -1) {
      sliceDomains[position] = i;
      position++;
    }
  }
 
  // Create new point data mpi group
  MPI_Group globalGroup, localGroup;
  MPI_Comm_group(mpiComm(), &globalGroup, AT_, "globalGroup");
  MPI_Group_incl(globalGroup, noSliceDomains, &sliceDomains[0], &localGroup, AT_);
 
  // Create new communicator and clean up
  MPI_Comm_create(mpiComm(), localGroup, &sliceSeries.m_mpiComm, AT_, "sliceSeries.m_mpiComm");
 
  MPI_Group_free(&globalGroup, AT_);
  MPI_Group_free(&localGroup, AT_);
 
  // Step 3: Exchange hilbertId data and determine offsets
 
  if(sliceSeries.m_noCellInSlice > 0) {
    // Get slice domain
    MInt sliceDomain = -1;
    MPI_Comm_rank(sliceSeries.m_mpiComm, &sliceDomain);
 
    // Determine total no of slice hilbertIds
    // Cells and elements are sorted by slice hilbertIds and also by domain. So if we have hilbertId
    // 2 on domain 3 and 5, then the order is determined by domaindId. This information is
    // distributed globally and the purpose of the next lines.
    MIntScratchSpace noLocalHilbertIds(noSliceDomains, AT_, "noLocalHilbertIds");
    noLocalHilbertIds[sliceDomain] = sliceSeries.m_noHilbertIds;
 
    MPI_Allgather(MPI_IN_PLACE, 1, MPI_INT, &noLocalHilbertIds[0], 1, MPI_INT, sliceSeries.m_mpiComm, AT_,
                  "MPI_IN_PLACE", "noLocalHilbertIds[0]");
 
    // Compute offset for recieve data
    MIntScratchSpace offsetsRecvData(noSliceDomains, AT_, "offsetsRecvData");
    offsetsRecvData[0] = 0;
    for(MInt i = 1; i < noSliceDomains; i++) {
      offsetsRecvData[i] = offsetsRecvData[i - 1] + noLocalHilbertIds[i - 1] * 3;
    }
    MInt noTotalRecvData = offsetsRecvData[noSliceDomains - 1] + noLocalHilbertIds[noSliceDomains - 1] * 3;
 
    // Create data array for sending
    MIntScratchSpace sendHilbertData(noLocalHilbertIds[sliceDomain] * 3, AT_, "sendHilbertData");
    // Send hilbertId, domainId and number of Nodes per HilbertId
    for(MInt i = 0; i < sliceSeries.m_noHilbertIds; i++) {
      sendHilbertData[i * 3] = hilbertInfo[i * 3];
      sendHilbertData[i * 3 + 1] = domainId();
      sendHilbertData[i * 3 + 2] = sliceSeries.m_elementNodesPerHilbertIds[i];
    }
 
    // Create data array for recieving
    MIntScratchSpace recvHilbertData(noTotalRecvData, AT_, "recvHilbertData");
    MIntScratchSpace noRecvData(noSliceDomains, AT_, "noRecvData");
    for(MInt i = 0; i < noSliceDomains; i++) {
      noRecvData[i] = noLocalHilbertIds[i] * 3;
    }
 
    // Exhange elementNoHilbertNodes
    MPI_Allgatherv(&sendHilbertData[0], sendHilbertData.size(), MPI_INT, &recvHilbertData[0], &noRecvData[0],
                   &offsetsRecvData[0], MPI_INT, sliceSeries.m_mpiComm, AT_, "sendHilbertData[0]",
                   "recvHilbertData[0]");
 
    // Create sliceGlobalHilbertInfo
    ScratchSpace<std::array<MInt, 3>> sliceGlobalHilbertInfo(noTotalRecvData / 3, AT_, "sliceGlobalHilbertInfo");
    for(MInt i = 0; i < noTotalRecvData / 3; i++) {
      // Contains number of cells sorted by hilbertId (first), domaindId (second) and number of
      // nodes per hilbertId
      sliceGlobalHilbertInfo[i][0] = recvHilbertData[i * 3];
      sliceGlobalHilbertInfo[i][1] = recvHilbertData[i * 3 + 1];
      sliceGlobalHilbertInfo[i][2] = recvHilbertData[i * 3 + 2];
    }
    // Sort sliceGlobalHilbertInfo by domainId...
    std::stable_sort(sliceGlobalHilbertInfo.begin(), sliceGlobalHilbertInfo.end(),
                     [](const std::array<MInt, 3>& a, const std::array<MInt, 3>& b) { return a[1] < b[1]; });
    // ... and then sort stable by hilbertId -> list ordered by hilbertId and then nested for each hilbertId by domainId
    std::stable_sort(sliceGlobalHilbertInfo.begin(), sliceGlobalHilbertInfo.end(),
                     [](const std::array<MInt, 3>& a, const std::array<MInt, 3>& b) { return a[0] < b[0]; });
 
    // Calulate elementNodeOffset by sliceHilbertIds for file writing
    MInt maxNoHilbertIds = *std::max_element(noLocalHilbertIds.begin(), noLocalHilbertIds.end());
    sliceSeries.m_elementNodesHilbertOffset.resize(maxNoHilbertIds);
    sliceSeries.m_cellIdsPerHilbertId.resize(maxNoHilbertIds);
    sliceSeries.m_cellIdsOffset.resize(maxNoHilbertIds);
    // The data in sliceGlobalHilbertInfo is sorted by their hilbertId on the first level. On the
    // second level by domainId. Example:
    //                             [hilbertId, domaindId, noNodes]
    // sliceGlobalHilbertInfo[0] = [0, 0, 20]
    // sliceGlobalHilbertInfo[1] = [0, 2, 10]
    // sliceGlobalHilbertInfo[2] = [1, 1, 6]
    //
    // To determine the global offset for the local data (for each local hilbertId) count the number
    // of nodes of all hilbertIds below the current one. If this hilbertId is distributed on
    // multiple domains, then count also the next number of nodes till the local domainid is reached
    // in sliceGlobalHilbertInfo. In the example above the offset for hilbertId 0 on domain 1 would
    // be 20.
    for(MInt i = 0, j = 0; i < noLocalHilbertIds[sliceDomain]; i++) {
      MInt offset = 0;
      j = 0;
      // Count until current hilbertId is reached
      while(sliceGlobalHilbertInfo[j][0] < hilbertInfo[i * 3]) {
        offset += sliceGlobalHilbertInfo[j][2];
        j++;
      }
      // Count to current domain
      while(sliceGlobalHilbertInfo[j][1] < domainId()) {
        offset += sliceGlobalHilbertInfo[j][2];
        j++;
      }
      // Set local node offset
      sliceSeries.m_elementNodesHilbertOffset[i] = offset;
      // save noCells with same hilbertId  for writing poly degree to the grid
      sliceSeries.m_cellIdsPerHilbertId[i] = hilbertInfo[i * 3 + 1];
      // save offset for writing poly degree to the grid
      sliceSeries.m_cellIdsOffset[i] = hilbertInfo[i * 3 + 2];
    }
  }
 
  // TODO labels:DG,DOC
  MBool optimizedSliceIo = true;
  optimizedSliceIo = Context::getBasicProperty<MBool>("optimizedSliceIo", AT_, &optimizedSliceIo);
 
  // Reorganize data before writing
  if(optimizedSliceIo && sliceSeries.m_noCellInSlice > 0) {
    // Create buffers to organize m_elementOffset
    std::vector<MInt> bufferElementIdsPerHilbertId;
    std::vector<MInt> bufferElementNodesPerHilbertIds;
    // Fill buffers
    std::copy(sliceSeries.m_elementIdsPerHilbertId.begin(),
              sliceSeries.m_elementIdsPerHilbertId.end(),
              back_inserter(bufferElementIdsPerHilbertId));
    std::copy(sliceSeries.m_elementNodesPerHilbertIds.begin(),
              sliceSeries.m_elementNodesPerHilbertIds.end(),
              back_inserter(bufferElementNodesPerHilbertIds));
 
    MInt noHilbertIdChunks = 0;
 
    // Reconstruct data for writing
    for(MInt h = 1, i = 0, offsetSum = 0; h <= (MInt)sliceSeries.m_elementNodesHilbertOffset.size(); h++) {
      // Reorganize data
      if(h < sliceSeries.m_noHilbertIds) {
        if((sliceSeries.m_elementNodesPerHilbertIds[i] + sliceSeries.m_elementNodesHilbertOffset[i])
           == sliceSeries.m_elementNodesHilbertOffset[h]) {
          sliceSeries.m_elementNodesPerHilbertIds[i] += sliceSeries.m_elementNodesPerHilbertIds[h];
          sliceSeries.m_elementIdsPerHilbertId[i] += sliceSeries.m_elementIdsPerHilbertId[h];
 
          sliceSeries.m_cellIdsPerHilbertId[i] += sliceSeries.m_cellIdsPerHilbertId[h];
 
          // Update element offsets of the elements belonging to the added hilbert id
          MInt startIdx = std::accumulate(&bufferElementIdsPerHilbertId[0], &bufferElementIdsPerHilbertId[h], 0);
          MInt endIdx = std::accumulate(&bufferElementIdsPerHilbertId[0], &bufferElementIdsPerHilbertId[h + 1], 0);
          offsetSum += bufferElementNodesPerHilbertIds[h - 1];
          for(MInt j = startIdx; j < endIdx; j++) {
            sliceSeries.m_elementOffset[j] += offsetSum;
          }
        } else {
          i++;
          sliceSeries.m_elementNodesPerHilbertIds[i] = sliceSeries.m_elementNodesPerHilbertIds[h];
          sliceSeries.m_elementNodesHilbertOffset[i] = sliceSeries.m_elementNodesHilbertOffset[h];
          sliceSeries.m_elementIdsPerHilbertId[i] = sliceSeries.m_elementIdsPerHilbertId[h];
 
          sliceSeries.m_cellIdsPerHilbertId[i] = sliceSeries.m_cellIdsPerHilbertId[h];
          sliceSeries.m_cellIdsOffset[i] = sliceSeries.m_cellIdsOffset[h];
 
          // Reset offsetSum in case of multiple concatenations
          offsetSum = 0;
        }
      }
 
      // Set remaining elements to zero at last iteration
      if(h == (MInt)sliceSeries.m_elementNodesHilbertOffset.size()) {
        i++;
        noHilbertIdChunks = i;
 
        std::fill(sliceSeries.m_elementNodesPerHilbertIds.begin() + i, sliceSeries.m_elementNodesPerHilbertIds.end(),
                  0);
        std::fill(sliceSeries.m_elementNodesHilbertOffset.begin() + i, sliceSeries.m_elementNodesHilbertOffset.end(),
                  0);
        std::fill(sliceSeries.m_elementIdsPerHilbertId.begin() + i, sliceSeries.m_elementIdsPerHilbertId.end(), 0);
 
        std::fill(sliceSeries.m_cellIdsPerHilbertId.begin() + i, sliceSeries.m_cellIdsPerHilbertId.end(), 0);
        std::fill(sliceSeries.m_cellIdsOffset.begin() + i, sliceSeries.m_cellIdsOffset.end(), 0);
      }
    }
 
    MInt maxNoHilbertIdChunks = -1;
    MPI_Allreduce(&noHilbertIdChunks, &maxNoHilbertIdChunks, 1, maia::type_traits<MInt>::mpiType(), MPI_MAX,
                  sliceSeries.m_mpiComm, AT_, "noHilbertIdChunks", "maxNoHilbertIdChunks");
    sliceSeries.m_elementNodesHilbertOffset.resize(maxNoHilbertIdChunks);
 
    sliceSeries.m_noHilbertIds = noHilbertIdChunks;
    sliceSeries.m_elementNodesPerHilbertIds.resize(noHilbertIdChunks);
    sliceSeries.m_elementIdsPerHilbertId.resize(noHilbertIdChunks);
    sliceSeries.m_cellIdsPerHilbertId.resize(noHilbertIdChunks);
    sliceSeries.m_cellIdsOffset.resize(noHilbertIdChunks);
  }
 
  sliceSeries.m_isInit = true;
  m_log << "done!" << std::endl;
}

mpiComm()

template<MInt nDim, class SysEqn >

MPI_Comm DgSlices< nDim, SysEqn >::mpiComm ( ) const

inlineprivate

Definition at line 36 of file dgcartesianslices.h.

{ return m_solver.mpiComm(); }

save() [1/2]

template<MInt nDim, class SysEqn >

void DgSlices< nDim, SysEqn >::save

(

const MBool

isFinalTimeStep

)

Author

Marcus Wiens (marcus) marcu.nosp@m.s.wi.nosp@m.ens@r.nosp@m.wth-.nosp@m.aache.nosp@m.n.de

Date

20.06.2016

Definition at line 761 of file dgcartesianslices.h.

                                                             {
  TRACE();
  // Return if dg slice is deactivated
  if(!m_enabled) {
    return;
  }
 
  // Call save method for all slice objects
  for(MUint i = 0; i < m_sliceSeries.size(); i++) {
    save(m_sliceSeries[i], isFinalTimeStep);
  }
}

save() [2/2]

template<MInt nDim, class SysEqn >

void DgSlices< nDim, SysEqn >::save ( DgSliceSeries &  sliceSeries, const MBool  isFinalTimeStep  )

private

Author

Marcus Wiens (marcus) marcu.nosp@m.s.wi.nosp@m.ens@r.nosp@m.wth-.nosp@m.aache.nosp@m.n.de

Date

06.01.2016

Definition at line 783 of file dgcartesianslices.h.

                                                                                         {
  TRACE();
 
  // Ensure that initalization was done
  if(!sliceSeries.m_isInit) {
    TERMM(1, "DG Slice class was not properly initialized.");
  }
 
  // Leave when domain is not part of slice
  if(sliceSeries.m_localNoNodes < 1) {
    return;
  }
 
  // Check for write Interval
  if(!(m_solver.m_timeStep % sliceSeries.m_writeInterval == 0 || isFinalTimeStep)) {
    return;
  }
 
  // TODO labels:DG @ansgar_slice not required at the moment
  /* // Compute sampling variables if not done already for another time series in this time step */
  /* if(m_solverCalcSamplingVars) { */
  /*   solver().calcSamplingVariables(m_solverSamplingVarIds); */
  /*   m_solverCalcSamplingVars = false; */
  /* } */
 
  // Write all slice variables to the same file
  const MInt noVars = std::accumulate(m_noSliceVars.begin(), m_noSliceVars.end(), 0);
  const MUint noVarIds = m_noSliceVars.size();
 
  // Get state at slice coordinates
  MFloatTensor coordinatesTensor(&sliceSeries.m_coordinates[0],
                                 std::ceil((MFloat)sliceSeries.m_coordinates.size() / m_dim), m_dim);
 
  // Save states of nodes
  MFloatScratchSpace states(std::ceil((MFloat)sliceSeries.m_coordinates.size() / m_dim) * noVars, AT_, "states");
  MFloatTensor stateTensor(&states[0], std::ceil((MFloat)sliceSeries.m_coordinates.size() / m_dim), noVars);
  // TODO labels:DG,totest check performance of state calculation (timer)
  // Calculate states of all nodes
  if(sliceSeries.m_coordinates.size() > 1) {
    for(MUint e = 0, offset = 0; e < sliceSeries.m_elementNodes.size(); e++) {
      const MUint noNodes = sliceSeries.m_elementNodes[e];
      for(MUint n = 0; n < noNodes; n++) {
        MInt varOffset = 0;
        for(MUint v = 0; v < noVarIds; v++) {
          m_solver.calcSamplingVarAtPoint(&coordinatesTensor(offset + n, 0), sliceSeries.m_elementIds[e],
                                          m_sliceVarIds[v], &stateTensor(offset + n, varOffset), true);
          varOffset += m_noSliceVars[v];
        }
      }
      offset += noNodes;
    }
  }
 
  // Write File
  using namespace maia::parallel_io;
 
  // TODO labels:DG slice naming with slice axis and intercept? same as in postprocessing
  // Create output file name with slices axis and number
  // TODO labels:DG change bX (block X) into e.g. sX (solver X)? need to change other output files as well
  const MString s = (m_solver.grid().raw().treeb().noSolvers() > 1) ? "_b" + std::to_string(solverId()) : "";
  std::ostringstream fileName;
  fileName << m_solver.outputDir() << "slice" << s << "_" << sliceSeries.m_fileNo << "_" << std::setw(8)
           << std::setfill('0') << m_solver.m_timeStep << ParallelIo::fileExt();
  ParallelIo file(fileName.str(), PIO_REPLACE, sliceSeries.m_mpiComm);
 
  // Determine offset and global number of nodes
  ParallelIo::size_type nodesOffset, globalNoNodes, offset, totalCount;
  ParallelIo::calcOffset(sliceSeries.m_localNoNodes, &nodesOffset, &globalNoNodes, sliceSeries.m_mpiComm);
  ParallelIo::calcOffset(sliceSeries.m_noCellInSlice, &offset, &totalCount, sliceSeries.m_mpiComm);
 
  // Set Attributes
  // Grid file name and solver type
  file.setAttribute(sliceSeries.m_gridFileName + ParallelIo::fileExt(), "gridFile");
  file.setAttribute(solverId(), "solverId");
  file.setAttribute("DG", "solverType");
  file.setAttribute(m_solver.m_timeStep, "timeStep");
  file.setAttribute(m_solver.m_time, "time");
  // Save slice axis and intercept to identify the grid file
  file.setAttribute(sliceSeries.m_axis, "sliceAxis");
  file.setAttribute(sliceSeries.m_intercept, "sliceIntercept");
 
  // Define arrays in file
  // Polyomial degree
  file.defineArray(PIO_INT, "polyDegs", totalCount);
 
  // Get information about integration method and polynomial type
  MString dgIntegrationMethod = "DG_INTEGRATE_GAUSS";
  dgIntegrationMethod =
      Context::getSolverProperty<MString>("dgIntegrationMethod", solverId(), AT_, &dgIntegrationMethod);
  MString dgPolynomialType = "DG_POLY_LEGENDRE";
  dgPolynomialType = Context::getSolverProperty<MString>("dgPolynomialType", solverId(), AT_, &dgPolynomialType);
 
  // Add integration method & polynomial type to grid file as attributes
  file.setAttribute(dgIntegrationMethod, "dgIntegrationMethod", "polyDegs");
  file.setAttribute(dgPolynomialType, "dgPolynomialType", "polyDegs");
 
  // Set Variables
  for(MUint v = 0, totalNoVars = 0; v < noVarIds; v++) {
    for(MInt i = 0; i < m_noSliceVars[v]; i++) {
      const MString name = "variables" + std::to_string(totalNoVars);
      file.defineArray(PIO_FLOAT, name, globalNoNodes);
      file.setAttribute(m_sliceVarNames[v][i], "name", name);
      totalNoVars++;
    }
  }
 
  // Write poly degree to file (sorted by hilbert id)
  for(MInt h = 0, pOffset = 0; h < (MInt)sliceSeries.m_elementNodesHilbertOffset.size(); h++) {
    if(h < sliceSeries.m_noHilbertIds && sliceSeries.m_cellIdsPerHilbertId[h] > 0) {
      file.setOffset(sliceSeries.m_cellIdsPerHilbertId[h], sliceSeries.m_cellIdsOffset[h]);
      file.writeArray(&sliceSeries.m_polyDegs[0 + pOffset], "polyDegs");
      pOffset += sliceSeries.m_cellIdsPerHilbertId[h];
    } else {
      file.setOffset(0, 0);
      file.writeArray(&sliceSeries.m_polyDegs[0], "polyDegs");
    }
  }
 
  // Number of nodes for elements per HilbertId
  const MInt tmp = std::accumulate(sliceSeries.m_elementNodes.begin(), sliceSeries.m_elementNodes.end(), 0);
  // create statetensor variable for easier accessing
  MFloatTensor stateTensorFinal(&states[0], std::max(tmp, 1), noVars);
 
  // Write data
  for(MInt i = 0; i < noVars; i++) {
    for(MInt h = 0, elementOffset = 0, offsetElementNodes = 0; h < (MInt)sliceSeries.m_elementNodesHilbertOffset.size();
        h++) {
      if(h < sliceSeries.m_noHilbertIds && sliceSeries.m_elementNodesPerHilbertIds[h] > 0) {
        // Write data for every (chunk of continuous) hilbertId(s) on this domain
 
        // Set amount and offset for current hilbertId
        file.setOffset(sliceSeries.m_elementNodesPerHilbertIds[h], sliceSeries.m_elementNodesHilbertOffset[h]);
        // Create buffer to collect data in correct order for saving
        MFloatScratchSpace buffer(sliceSeries.m_elementNodesPerHilbertIds[h], AT_, "buffer");
        std::fill(buffer.begin(), buffer.end(), 0.0);
        const MInt noElementsSlice = sliceSeries.m_elementIdsPerHilbertId[h] + elementOffset;
 
        // Fill the buffer at positions where the element data belongs (reminder: elements only
        // exist on highest refinement level)
        for(MInt e = elementOffset; e < noElementsSlice; e++) {
          MFloat* const b = &buffer[sliceSeries.m_elementOffset[e]];
          for(MInt j = 0; j < sliceSeries.m_elementNodes[e]; j++) {
            b[j] = stateTensorFinal(offsetElementNodes + j, i);
          }
          offsetElementNodes += sliceSeries.m_elementNodes[e];
        }
        elementOffset = noElementsSlice;
 
        const MString name = "variables" + std::to_string(i);
        // write data
        file.writeArray(&buffer[0], name);
      } else {
        // Dummy calls for data writing if this domain is finished, but other domains not
        file.setOffset(0, 0);
        const MFloat tmpBuf = 0;
        const MString name = "variables" + std::to_string(i);
        file.writeArray(&tmpBuf, name);
      }
    }
  }
}

setProperties()

template<MInt nDim, class SysEqn >

void DgSlices< nDim, SysEqn >::setProperties

Author

Marcus Wiens (marcus) marcu.nosp@m.s.wi.nosp@m.ens@r.nosp@m.wth-.nosp@m.aache.nosp@m.n.de

Date

06.01.2016

Definition at line 138 of file dgcartesianslices.h.

solverId()

template<MInt nDim, class SysEqn >

MInt DgSlices< nDim, SysEqn >::solverId ( ) const

inlineprivate

Definition at line 37 of file dgcartesianslices.h.

{ return m_solver.m_solverId; }

updateGridFile()

template<MInt nDim, class SysEqn >

void DgSlices< nDim, SysEqn >::updateGridFile ( DgSliceSeries &  sliceSeries, const MString &  prefix = ""  )

private

Member Data Documentation

m_dim

template<MInt nDim, class SysEqn >

constexpr const MInt DgSlices< nDim, SysEqn >::m_dim = 3

staticconstexprprivate

Definition at line 46 of file dgcartesianslices.h.

m_enabled

template<MInt nDim, class SysEqn >

MBool DgSlices< nDim, SysEqn >::m_enabled = false

private

Definition at line 51 of file dgcartesianslices.h.

m_noSliceVars

template<MInt nDim, class SysEqn >

std::vector<MInt> DgSlices< nDim, SysEqn >::m_noSliceVars {}

private

Definition at line 58 of file dgcartesianslices.h.

m_sliceSeries

template<MInt nDim, class SysEqn >

std::vector<DgSliceSeries> DgSlices< nDim, SysEqn >::m_sliceSeries {}

private

Definition at line 53 of file dgcartesianslices.h.

m_sliceVarIds

template<MInt nDim, class SysEqn >

std::vector<MInt> DgSlices< nDim, SysEqn >::m_sliceVarIds {}

private

Definition at line 56 of file dgcartesianslices.h.

m_sliceVarNames

template<MInt nDim, class SysEqn >

std::vector<std::vector<MString> > DgSlices< nDim, SysEqn >::m_sliceVarNames {}

private

Definition at line 60 of file dgcartesianslices.h.

m_solver

template<MInt nDim, class SysEqn >

DgCartesianSolver<nDim, SysEqn>& DgSlices< nDim, SysEqn >::m_solver

private

Definition at line 49 of file dgcartesianslices.h.

---

The documentation for this class was generated from the following file:

• /home/gitlab-runner/scratch/builds/oxpnswJ6/1/aia/m-AIA/m-AIA/src/DG/dgcartesianslices.h