Cmfd refactor to sync with ClosedMOC/3D-MOC

src/CPUSolver.cpp

@@ -641,7 +641,7 @@ void CPUSolver::computeKeff() {
 
  /* Reduce new fission rates across FSRs */
  fission = pairwise_sum<FP_PRECISION>(FSR_rates, _num_FSRs);
- 
+
  _k_eff *= fission;
 
  delete [] FSR_rates;
@@ -775,7 +775,7 @@ void CPUSolver::tallyScalarFlux(segment* curr_segment, int azim_index,
 
 /**
  * @brief Tallies the current contribution from this segment across the
- * the appropriate CMFD mesh cell surface or corner.
+ * the appropriate CMFD mesh cell surface.
  * @param curr_segment a pointer to the Track segment of interest
  * @param azim_index the azimuthal index for this segmenbt
  * @param track_flux a pointer to the Track's angular flux
@@ -784,7 +784,7 @@ void CPUSolver::tallyScalarFlux(segment* curr_segment, int azim_index,
 void CPUSolver::tallyCurrent(segment* curr_segment, int azim_index,
  FP_PRECISION* track_flux, bool fwd) {
 
- /* Tally surface or corner currents if CMFD is in use */
+ /* Tally surface currents if CMFD is in use */
  if (_cmfd != NULL && _cmfd->isFluxUpdateOn())
  _cmfd->tallyCurrent(curr_segment, track_flux,
  &_polar_weights(azim_index,0), fwd);

src/CPUSolver.h

@@ -55,8 +55,7 @@ class CPUSolver : public Solver {
  FP_PRECISION* track_flux, FP_PRECISION* fsr_flux);
 
  /**
- * @brief Computes the contribution to surface or corner current from a Track
- * segment.
+ * @brief Computes the contribution to surface current from a segment.
  * @param curr_segment a pointer to the Track segment of interest
  * @param azim_index a pointer to the azimuthal angle index for this segment
  * @param track_flux a pointer to the Track's angular flux

src/Cmfd.h

@@ -137,9 +137,6 @@ class Cmfd {
  /** Vector of surface currents for each CMFD cell */
  Vector* _surface_currents;
 
- /** Vector of corner currents for each CMFD cell */
- Vector* _corner_currents;
-
  /** Vector of vectors of FSRs containing in each cell */
  std::vector< std::vector<int> > _cell_fsrs;
 
@@ -159,14 +156,18 @@ class Cmfd {
  std::map<int, std::vector< std::pair<int, FP_PRECISION> > >
  _k_nearest_stencils;
 
+ /** OpenMP mutual exclusion locks for atomic CMFD cell operations */
+ omp_lock_t* _cell_locks;
+
  /* Private worker functions */
  FP_PRECISION computeLarsensEDCFactor(FP_PRECISION dif_coef,
  FP_PRECISION delta);
  void constructMatrices(int moc_iteration);
  void collapseXS();
  void updateMOCFlux();
  void rescaleFlux();
- void splitCorners();
+ void splitEdgeCurrents();
+ void getEdgeSplitSurfaces(int cell, int edge, std::vector<int>* surfaces);
  void initializeMaterials();
  void initializeCurrents();
  void generateKNearestStencils();
@@ -198,7 +199,6 @@ class Cmfd {
  void initializeLattice(Point* offset);
  int findCmfdCell(LocalCoords* coords);
  int findCmfdSurface(int cell_id, LocalCoords* coords);
- int findCmfdCorner(int cell_id, LocalCoords* coords);
  void addFSRToCell(int cell_id, int fsr_id);
  void zeroCurrents();
  void tallyCurrent(segment* curr_segment, FP_PRECISION* track_flux,

src/Geometry.cpp

@@ -846,8 +846,6 @@ void Geometry::segmentize(Track* track) {
  new_segment->_cmfd_surface_fwd = _cmfd->findCmfdSurface(cmfd_cell, &end);
  new_segment->_cmfd_surface_bwd =
  _cmfd->findCmfdSurface(cmfd_cell, &start);
- new_segment->_cmfd_corner_fwd = _cmfd->findCmfdCorner(cmfd_cell, &end);
- new_segment->_cmfd_corner_bwd = _cmfd->findCmfdCorner(cmfd_cell, &start);
 
  /* Re-nudge segments from surface */
  start.adjustCoords(TINY_MOVE);
@@ -1070,10 +1068,9 @@ void Geometry::initializeCmfd() {
 
  /* Initialize the CMFD lattice */
  Point offset;
- double offset_x = getMinX() + getWidthX()/2.0;
- double offset_y = getMinY() + getWidthY()/2.0;
- offset.setX(offset_x);
- offset.setY(offset_y);
+ offset.setX(getMinX() + getWidthX()/2.0);
+ offset.setY(getMinY() + getWidthY()/2.0);
+
  _cmfd->initializeLattice(&offset);
 }
 

src/Material.cpp

@@ -420,7 +420,7 @@ void Material::setVolume(double volume) {
 
 /**
  * @brief Increment the volume/area of the Material by some amount.
- * @details This routine is called by the TrackGenerator during track 
+ * @details This routine is called by the TrackGenerator during track
  * generation and segmentation.
  * @param volume the amount to increment the current volume by
  */
@@ -440,7 +440,7 @@ void Material::setNumInstances(int num_instances) {
 
 /**
  * @brief Increment the number of instances of this Material.
- * @details This routine is called by the TrackGenerator during track 
+ * @details This routine is called by the TrackGenerator during track
  * generation and segmentation.
  */
 void Material::incrementNumInstances() {

src/Matrix.cpp

@@ -11,11 +11,12 @@
  * outside. Locks are used to make the matrix thread-safe against
  * concurrent writes the same value. One lock locks out multiple rows of
  * the matrix at a time reprsenting multiple groups in the same cell.
+ * @param cell_locks Omp locks for atomic cell operations
  * @param num_x The number of cells in the x direction.
  * @param num_y The number of cells in the y direction.
  * @param num_groups The number of energy groups in each cell.
  */
-Matrix::Matrix(int num_x, int num_y, int num_groups) {
+Matrix::Matrix(omp_lock_t* cell_locks, int num_x, int num_y, int num_groups) {
 
  setNumX(num_x);
  setNumY(num_y);
@@ -32,13 +33,12 @@ Matrix::Matrix(int num_x, int num_y, int num_groups) {
  _DIAG = NULL;
  _modified = true;
 
- /* Allocate memory for OpenMP locks for each Matrix cell */
- _cell_locks = new omp_lock_t[_num_x*_num_y];
+ /* Set OpenMP locks for each Matrix cell */
+ if (cell_locks == NULL)
+ log_printf(ERROR, "Unable to create a Matrix without an array of cell "
+ "locks");
 
- /* Loop over all Matrix cells to initialize OpenMP locks */
- #pragma omp parallel for schedule(guided)
- for (int r=0; r < _num_x*_num_y; r++)
- omp_init_lock(&_cell_locks[r]);
+ _cell_locks = cell_locks;
 }
 
 
@@ -63,9 +63,6 @@ Matrix::~Matrix() {
  for (int i=0; i < _num_rows; i++)
  _LIL[i].clear();
  _LIL.clear();
-
- if (_cell_locks != NULL)
- delete [] _cell_locks;
 }
 
 
@@ -426,7 +423,7 @@ void Matrix::setNumGroups(int num_groups) {
  */
 void Matrix::transpose() {
 
- Matrix temp(_num_x, _num_y, _num_groups);
+ Matrix temp(_cell_locks, _num_x, _num_y, _num_groups);
  convertToCSR();
  int col, cell_to, cell_from, group_to, group_from;
  FP_PRECISION val;
@@ -463,3 +460,12 @@ void Matrix::transpose() {
  }
  }
 }
+
+
+/**
+ * @brief Return the array of cell locks for atomic cell operations.
+ * @return an array of cell locks
+ */
+omp_lock_t* Matrix::getCellLocks() {
+ return _cell_locks;
+}

src/Matrix.h

@@ -50,7 +50,7 @@ class Matrix {
  void setNumGroups(int num_groups);
 
 public:
- Matrix(int num_x=1, int num_y=1, int num_groups=1);
+ Matrix(omp_lock_t* cell_locks, int num_x=1, int num_y=1, int num_groups=1);
  virtual ~Matrix();
 
  /* Worker functions */
@@ -72,6 +72,7 @@ class Matrix {
  int getNumGroups();
  int getNumRows();
  int getNNZ();
+ omp_lock_t* getCellLocks();
 
  /* Setter functions */
  void setValue(int cell_from, int group_from, int cell_to, int group_to,

src/Solver.cpp

@@ -596,7 +596,7 @@ void Solver::initializeFSRs() {
  _FSR_volumes = _track_generator->getFSRVolumes();
 
  /* Generate the FSR centroids */
- _track_generator->generateFSRCentroids();
+ _track_generator->generateFSRCentroids(_FSR_volumes);
 
  /* Attach the correct materials to each track segment */
  _track_generator->initializeSegments();
@@ -834,7 +834,7 @@ void Solver::computeFlux(int max_iters, solverMode mode,
  _k_eff = 1.;
 
  _num_iterations = 0;
- FP_PRECISION residual;
+ FP_PRECISION residual = 0.;
 
  /* Initialize data structures */
  initializeFSRs();
@@ -944,7 +944,7 @@ void Solver::computeSource(int max_iters, solverMode mode,
  _k_eff = k_eff;
 
  _num_iterations = 0;
- FP_PRECISION residual;
+ FP_PRECISION residual = 0.;
 
  /* Initialize data structures */
  initializeFSRs();
@@ -1026,7 +1026,7 @@ void Solver::computeEigenvalue(int max_iters, solverMode mode,
  _timer->startTimer();
 
  _num_iterations = 0;
- FP_PRECISION residual;
+ FP_PRECISION residual = 0.;
 
  /* An initial guess for the eigenvalue */
  _k_eff = 1.0;

src/Track.h

@@ -41,18 +41,10 @@ struct segment {
  /** The ID for the mesh surface crossed by the Track start point */
  int _cmfd_surface_bwd;
 
- /** The ID for the mesh corner crossed by the Track end point */
- int _cmfd_corner_fwd;
-
- /** The ID for the mesh corner crossed by the Track start point */
- int _cmfd_corner_bwd;
-
  /** Constructor initializes CMFD surfaces */
  segment() {
  _cmfd_surface_fwd = -1;
  _cmfd_surface_bwd = -1;
- _cmfd_corner_fwd = -1;
- _cmfd_corner_bwd = -1;
  }
 };
 

src/TrackGenerator.cpp

@@ -1339,8 +1339,6 @@ void TrackGenerator::dumpTracksToFile() {
  int region_id;
  int cmfd_surface_fwd;
  int cmfd_surface_bwd;
- int cmfd_corner_fwd;
- int cmfd_corner_bwd;
 
  /* Loop over all Tracks */
  for (int i=0; i < _num_azim; i++) {
@@ -1387,12 +1385,8 @@ void TrackGenerator::dumpTracksToFile() {
  if (cmfd != NULL) {
  cmfd_surface_fwd = curr_segment->_cmfd_surface_fwd;
  cmfd_surface_bwd = curr_segment->_cmfd_surface_bwd;
- cmfd_corner_fwd = curr_segment->_cmfd_corner_fwd;
- cmfd_corner_bwd = curr_segment->_cmfd_corner_bwd;
  fwrite(&cmfd_surface_fwd, sizeof(int), 1, out);
  fwrite(&cmfd_surface_bwd, sizeof(int), 1, out);
- fwrite(&cmfd_corner_fwd, sizeof(int), 1, out);
- fwrite(&cmfd_corner_bwd, sizeof(int), 1, out);
  }
  }
  }
@@ -1555,8 +1549,6 @@ bool TrackGenerator::readTracksFromFile() {
 
  int cmfd_surface_fwd;
  int cmfd_surface_bwd;
- int cmfd_corner_fwd;
- int cmfd_corner_bwd;
  segment curr_segment;
 
  std::map<int, Material*> materials = _geometry->getAllMaterials();
@@ -1601,12 +1593,8 @@ bool TrackGenerator::readTracksFromFile() {
  if (cmfd != NULL) {
  ret = fread(&cmfd_surface_fwd, sizeof(int), 1, in);
  ret = fread(&cmfd_surface_bwd, sizeof(int), 1, in);
- ret = fread(&cmfd_corner_fwd, sizeof(int), 1, in);
- ret = fread(&cmfd_corner_bwd, sizeof(int), 1, in);
  curr_segment._cmfd_surface_fwd = cmfd_surface_fwd;
  curr_segment._cmfd_surface_bwd = cmfd_surface_bwd;
- curr_segment._cmfd_corner_fwd = cmfd_corner_fwd;
- curr_segment._cmfd_corner_bwd = cmfd_corner_bwd;
  }
 
  /* Add this segment to the Track */
@@ -1772,14 +1760,12 @@ void TrackGenerator::correctFSRVolume(int fsr_id, FP_PRECISION fsr_volume) {
  * FSR by the segment's length and azimuthal weight. The numerical
  * centroid fomula can be found in R. Ferrer et. al. "Linear Source
  * Approximation in CASMO 5", PHYSOR 2012.
+ * @param FSR_volumes An array of FSR volumes.
  */
-void TrackGenerator::generateFSRCentroids() {
+void TrackGenerator::generateFSRCentroids(FP_PRECISION* FSR_volumes) {
 
  int num_FSRs = _geometry->getNumFSRs();
 
- /* Get FSR Volumes */
- FP_PRECISION* FSR_volumes = getFSRVolumes();
-
  /* Create array of centroids and initialize to origin */
  Point** centroids = new Point*[num_FSRs];
  for (int r=0; r < num_FSRs; r++) {
@@ -1821,7 +1807,6 @@ void TrackGenerator::generateFSRCentroids() {
  _geometry->setFSRCentroid(r, centroids[r]);
 
  /* Delete temporary array of FSR volumes and centroids */
- delete [] FSR_volumes;
  delete [] centroids;
 }
 
@@ -1849,7 +1834,6 @@ void TrackGenerator::splitSegments(FP_PRECISION max_optical_length) {
  FP_PRECISION* sigma_t;
  int num_groups;
  int cmfd_surface_fwd, cmfd_surface_bwd;
- int cmfd_corner_fwd, cmfd_corner_bwd;
 
  /* Iterate over all Tracks */
  for (int i=0; i < _num_azim; i++) {
@@ -1863,8 +1847,6 @@ void TrackGenerator::splitSegments(FP_PRECISION max_optical_length) {
  fsr_id = curr_segment->_region_id;
  cmfd_surface_fwd = curr_segment->_cmfd_surface_fwd;
  cmfd_surface_bwd = curr_segment->_cmfd_surface_bwd;
- cmfd_corner_fwd = curr_segment->_cmfd_corner_fwd;
- cmfd_corner_bwd = curr_segment->_cmfd_corner_bwd;
 
  /* Compute number of segments to split this segment into */
  min_num_cuts = 1;
@@ -1891,15 +1873,11 @@ void TrackGenerator::splitSegments(FP_PRECISION max_optical_length) {
  new_segment->_region_id = fsr_id;
 
  /* Assign CMFD surface boundaries */
- if (k == 0) {
+ if (k == 0)
  new_segment->_cmfd_surface_bwd = cmfd_surface_bwd;
- new_segment->_cmfd_corner_bwd = cmfd_corner_bwd;
- }
 
- if (k == min_num_cuts-1) {
+ if (k == min_num_cuts-1)
  new_segment->_cmfd_surface_fwd = cmfd_surface_fwd;
- new_segment->_cmfd_corner_fwd = cmfd_corner_fwd;
- }
 
  /* Insert the new segment to the Track */
  _tracks[i][j].insertSegment(s+k+1, new_segment);

src/TrackGenerator.h

@@ -141,7 +141,7 @@ class TrackGenerator {
  void retrieveSegmentCoords(double* coords, int num_segments);
  void generateTracks(bool neighbor_cells=false);
  void correctFSRVolume(int fsr_id, FP_PRECISION fsr_volume);
- void generateFSRCentroids();
+ void generateFSRCentroids(FP_PRECISION* FSR_volumes);
  void splitSegments(FP_PRECISION max_optical_length);
  void initializeSegments();
 };