23084721708

Committed 14 Mar 2026 08:56AM UTC coverage: 81.599% (-0.5%) from 82.058%

Build # 23084721708

Build Type

Pull #2693

github

Committed by

web-flow

Commit Message

Merge 0ed23ee59 into bc9c31e0f

Pull Request Pull Request #2693: Add reactivity control to coupled transport-depletion analyses

Run Details

17575 of 25275 branches covered (69.54%)

Branch coverage included in aggregate %.

74 of 85 new or added lines in 4 files covered. (87.06%)

3755 existing lines in 99 files now uncovered.

58074 of 67433 relevant lines covered (86.12%)

47252067.37 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

91.64

/src/cmfd_solver.cpp

#include "openmc/cmfd_solver.h"

#include <cmath>

#ifdef _OPENMP
#include <omp.h>
#endif
#include "openmc/tensor.h"

#include "openmc/bank.h"
#include "openmc/capi.h"
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/mesh.h"
#include "openmc/message_passing.h"
#include "openmc/tallies/filter_energy.h"
#include "openmc/tallies/filter_mesh.h"
#include "openmc/tallies/tally.h"
#include "openmc/vector.h"

namespace openmc {

namespace cmfd {

//==============================================================================
// Global variables
//==============================================================================

vector<int> indptr;

vector<int> indices;

int dim;

double spectral;

int nx, ny, nz, ng;

tensor::Tensor<int> indexmap;

int use_all_threads;

StructuredMesh* mesh;

vector<double> egrid;

double norm;

} // namespace cmfd

//==============================================================================
// GET_CMFD_ENERGY_BIN returns the energy bin for a source site energy
//==============================================================================

int get_cmfd_energy_bin(const double E)
{
  // Check if energy is out of grid bounds
  if (E < cmfd::egrid[0]) {
    // throw warning message
    warning("Detected source point below energy grid");
    return 0;
  } else if (E >= cmfd::egrid[cmfd::ng]) {
    // throw warning message
    warning("Detected source point above energy grid");
    return cmfd::ng - 1;
  } else {
    // Iterate through energy grid to find matching bin
    for (int g = 0; g < cmfd::ng; g++) {
      if (E >= cmfd::egrid[g] && E < cmfd::egrid[g + 1]) {
        return g;
      }
    }
  }
  // Return -1 by default
  return -1;
}

//==============================================================================
// COUNT_BANK_SITES bins fission sites according to CMFD mesh and energy
//==============================================================================

tensor::Tensor<double> count_bank_sites(
  tensor::Tensor<int>& bins, bool* outside)
{
  // Determine shape of array for counts
  std::size_t cnt_size = cmfd::nx * cmfd::ny * cmfd::nz * cmfd::ng;

  // Create array of zeros
  tensor::Tensor<double> cnt = tensor::zeros<double>({cnt_size});
  bool outside_ = false;

  auto bank_size = simulation::source_bank.size();
  for (int i = 0; i < bank_size; i++) {
    const auto& site = simulation::source_bank[i];

    // determine scoring bin for CMFD mesh
    int mesh_bin = cmfd::mesh->get_bin(site.r);

    // if outside mesh, skip particle
    if (mesh_bin < 0) {
      outside_ = true;
      continue;
    }

    // determine scoring bin for CMFD energy
    int energy_bin = get_cmfd_energy_bin(site.E);

    // add to appropriate bin
    cnt(mesh_bin * cmfd::ng + energy_bin) += site.wgt;

    // store bin index which is used again when updating weights
    bins[i] = mesh_bin * cmfd::ng + energy_bin;
  }

  int total = cnt.size();
  tensor::Tensor<double> counts = tensor::zeros<double>({cnt_size});

#ifdef OPENMC_MPI
  // collect values from all processors
  MPI_Reduce(
    cnt.data(), counts.data(), total, MPI_DOUBLE, MPI_SUM, 0, mpi::intracomm);

  // Check if there were sites outside the mesh for any processor
  MPI_Reduce(&outside_, outside, 1, MPI_C_BOOL, MPI_LOR, 0, mpi::intracomm);

#else
  std::copy(cnt.data(), cnt.data() + total, counts.data());
  *outside = outside_;
#endif

  return counts;
}

//==============================================================================
// OPENMC_CMFD_REWEIGHT performs reweighting of particles in source bank
//==============================================================================

extern "C" void openmc_cmfd_reweight(
  const bool feedback, const double* cmfd_src)
{
  // Get size of source bank and cmfd_src
  auto bank_size = simulation::source_bank.size();
  std::size_t src_size = cmfd::nx * cmfd::ny * cmfd::nz * cmfd::ng;

  // count bank sites for CMFD mesh, store bins in bank_bins for reweighting
  tensor::Tensor<int> bank_bins = tensor::zeros<int>({bank_size});
  bool sites_outside;
  tensor::Tensor<double> sourcecounts =
    count_bank_sites(bank_bins, &sites_outside);

  // Compute CMFD weightfactors
  tensor::Tensor<double> weightfactors = tensor::ones<double>({src_size});
  if (mpi::master) {
    if (sites_outside) {
      fatal_error("Source sites outside of the CMFD mesh");
    }

    double norm = sourcecounts.sum() / cmfd::norm;
    for (int i = 0; i < src_size; i++) {
      if (sourcecounts[i] > 0 && cmfd_src[i] > 0) {
        weightfactors[i] = cmfd_src[i] * norm / sourcecounts[i];
      }
    }
  }

  if (!feedback)
    return;

#ifdef OPENMC_MPI
  // Send weightfactors to all processors
  MPI_Bcast(weightfactors.data(), src_size, MPI_DOUBLE, 0, mpi::intracomm);
#endif

  // Iterate through fission bank and update particle weights
  for (int64_t i = 0; i < bank_size; i++) {
    auto& site = simulation::source_bank[i];
    site.wgt *= weightfactors(bank_bins(i));
  }
}

//==============================================================================
// OPENMC_INITIALIZE_MESH_EGRID sets the mesh and energy grid for CMFD reweight
//==============================================================================

extern "C" void openmc_initialize_mesh_egrid(
  const int meshtally_id, const int* cmfd_indices, const double norm)
{
  // Make sure all CMFD memory is freed
  free_memory_cmfd();

  // Set CMFD indices
  cmfd::nx = cmfd_indices[0];
  cmfd::ny = cmfd_indices[1];
  cmfd::nz = cmfd_indices[2];
  cmfd::ng = cmfd_indices[3];

  // Set CMFD reweight properties
  cmfd::norm = norm;

  // Find index corresponding to tally id
  int32_t tally_index;
  openmc_get_tally_index(meshtally_id, &tally_index);

  // Get filters assocaited with tally
  const auto& tally_filters = model::tallies[tally_index]->filters();

  // Get mesh filter index
  auto meshfilter_index = tally_filters[0];

  // Store energy filter index if defined, otherwise set to -1
  auto energy_index = (tally_filters.size() == 2) ? tally_filters[1] : -1;

  // Get mesh index from mesh filter index
  int32_t mesh_index;
  openmc_mesh_filter_get_mesh(meshfilter_index, &mesh_index);

  // Get mesh from mesh index
  cmfd::mesh = dynamic_cast<StructuredMesh*>(model::meshes[mesh_index].get());

  // Get energy bins from energy index, otherwise use default
  if (energy_index != -1) {
    auto efilt_base = model::tally_filters[energy_index].get();
    auto* efilt = dynamic_cast<EnergyFilter*>(efilt_base);
    cmfd::egrid = efilt->bins();
  } else {
    cmfd::egrid = {0.0, INFTY};
  }
}

//==============================================================================
// MATRIX_TO_INDICES converts a matrix index to spatial and group
// indices
//==============================================================================

void matrix_to_indices(int irow, int& g, int& i, int& j, int& k)
{
  g = irow % cmfd::ng;
  i = cmfd::indexmap(irow / cmfd::ng, 0);
  j = cmfd::indexmap(irow / cmfd::ng, 1);
  k = cmfd::indexmap(irow / cmfd::ng, 2);
}

//==============================================================================
// GET_DIAGONAL_INDEX returns the index in CSR index array corresponding to
// the diagonal element of a specified row
//==============================================================================

int get_diagonal_index(int row)
{
  for (int j = cmfd::indptr[row]; j < cmfd::indptr[row + 1]; j++) {
    if (cmfd::indices[j] == row)
      return j;
  }

  // Return -1 if not found
  return -1;
}

//==============================================================================
// SET_INDEXMAP sets the elements of indexmap based on input coremap
//==============================================================================

void set_indexmap(const int* coremap)
{
  for (int z = 0; z < cmfd::nz; z++) {
    for (int y = 0; y < cmfd::ny; y++) {
      for (int x = 0; x < cmfd::nx; x++) {
        int idx = (z * cmfd::ny * cmfd::nx) + (y * cmfd::nx) + x;
        if (coremap[idx] != CMFD_NOACCEL) {
          int counter = coremap[idx];
          cmfd::indexmap(counter, 0) = x;
          cmfd::indexmap(counter, 1) = y;
          cmfd::indexmap(counter, 2) = z;
        }
      }
    }
  }
}

//==============================================================================
// CMFD_LINSOLVER_1G solves a one group CMFD linear system
//==============================================================================

int cmfd_linsolver_1g(
  const double* A_data, const double* b, double* x, double tol)
{
  // Set overrelaxation parameter
  double w = 1.0;

  // Perform Gauss-Seidel iterations
  for (int igs = 1; igs <= 10000; igs++) {
    double err = 0.0;

    // Copy over x vector
    vector<double> tmpx {x, x + cmfd::dim};

    // Perform red/black Gauss-Seidel iterations
    for (int irb = 0; irb < 2; irb++) {

// Loop around matrix rows
#pragma omp parallel for reduction(+ : err) if (cmfd::use_all_threads)
      for (int irow = 0; irow < cmfd::dim; irow++) {
        int g, i, j, k;
        matrix_to_indices(irow, g, i, j, k);

        // Filter out black cells
        if ((i + j + k) % 2 != irb)
          continue;

        // Get index of diagonal for current row
        int didx = get_diagonal_index(irow);

        // Perform temporary sums, first do left of diag, then right of diag
        double tmp1 = 0.0;
        for (int icol = cmfd::indptr[irow]; icol < didx; icol++)
          tmp1 += A_data[icol] * x[cmfd::indices[icol]];
        for (int icol = didx + 1; icol < cmfd::indptr[irow + 1]; icol++)
          tmp1 += A_data[icol] * x[cmfd::indices[icol]];

        // Solve for new x
        double x1 = (b[irow] - tmp1) / A_data[didx];

        // Perform overrelaxation
        x[irow] = (1.0 - w) * x[irow] + w * x1;

        // Compute residual and update error
        double res = (tmpx[irow] - x[irow]) / tmpx[irow];
        err += res * res;
      }
    }

    // Check convergence
    err = std::sqrt(err / cmfd::dim);
    if (err < tol)
      return igs;

    // Calculate new overrelaxation parameter
    w = 1.0 / (1.0 - 0.25 * cmfd::spectral * w);
  }

  // Throw error, as max iterations met
  fatal_error("Maximum Gauss-Seidel iterations encountered.");

  // Return -1 by default, although error thrown before reaching this point
  return -1;
}

//==============================================================================
// CMFD_LINSOLVER_2G solves a two group CMFD linear system
//==============================================================================

int cmfd_linsolver_2g(
  const double* A_data, const double* b, double* x, double tol)
{
  // Set overrelaxation parameter
  double w = 1.0;

  // Perform Gauss-Seidel iterations
  for (int igs = 1; igs <= 10000; igs++) {
    double err = 0.0;

    // Copy over x vector
    vector<double> tmpx {x, x + cmfd::dim};

    // Perform red/black Gauss-Seidel iterations
    for (int irb = 0; irb < 2; irb++) {

// Loop around matrix rows
#pragma omp parallel for reduction(+ : err) if (cmfd::use_all_threads)
      for (int irow = 0; irow < cmfd::dim; irow += 2) {
        int g, i, j, k;
        matrix_to_indices(irow, g, i, j, k);

        // Filter out black cells
        if ((i + j + k) % 2 != irb)
          continue;

        // Get index of diagonals for current row and next row
        int d1idx = get_diagonal_index(irow);
        int d2idx = get_diagonal_index(irow + 1);

        // Get block diagonal
        double m11 = A_data[d1idx]; // group 1 diagonal
        double m12 =
          A_data[d1idx + 1]; // group 1 right of diagonal (sorted by col)
        double m21 =
          A_data[d2idx - 1];        // group 2 left of diagonal (sorted by col)
        double m22 = A_data[d2idx]; // group 2 diagonal

        // Analytically invert the diagonal
        double dm = m11 * m22 - m12 * m21;
        double d11 = m22 / dm;
        double d12 = -m12 / dm;
        double d21 = -m21 / dm;
        double d22 = m11 / dm;

        // Perform temporary sums, first do left of diag, then right of diag
        double tmp1 = 0.0;
        double tmp2 = 0.0;
        for (int icol = cmfd::indptr[irow]; icol < d1idx; icol++)
          tmp1 += A_data[icol] * x[cmfd::indices[icol]];
        for (int icol = cmfd::indptr[irow + 1]; icol < d2idx - 1; icol++)
          tmp2 += A_data[icol] * x[cmfd::indices[icol]];
        for (int icol = d1idx + 2; icol < cmfd::indptr[irow + 1]; icol++)
          tmp1 += A_data[icol] * x[cmfd::indices[icol]];
        for (int icol = d2idx + 1; icol < cmfd::indptr[irow + 2]; icol++)
          tmp2 += A_data[icol] * x[cmfd::indices[icol]];

        // Adjust with RHS vector
        tmp1 = b[irow] - tmp1;
        tmp2 = b[irow + 1] - tmp2;

        // Solve for new x
        double x1 = d11 * tmp1 + d12 * tmp2;
        double x2 = d21 * tmp1 + d22 * tmp2;

        // Perform overrelaxation
        x[irow] = (1.0 - w) * x[irow] + w * x1;
        x[irow + 1] = (1.0 - w) * x[irow + 1] + w * x2;

        // Compute residual and update error
        double res = (tmpx[irow] - x[irow]) / tmpx[irow];
        err += res * res;
      }
    }

    // Check convergence
    err = std::sqrt(err / cmfd::dim);
    if (err < tol)
      return igs;

    // Calculate new overrelaxation parameter
    w = 1.0 / (1.0 - 0.25 * cmfd::spectral * w);
  }

  // Throw error, as max iterations met
  fatal_error("Maximum Gauss-Seidel iterations encountered.");

  // Return -1 by default, although error thrown before reaching this point
  return -1;
}

//==============================================================================
// CMFD_LINSOLVER_NG solves a general CMFD linear system
//==============================================================================

int cmfd_linsolver_ng(
  const double* A_data, const double* b, double* x, double tol)
{
  // Set overrelaxation parameter
  double w = 1.0;

  // Perform Gauss-Seidel iterations
  for (int igs = 1; igs <= 10000; igs++) {
    double err = 0.0;

    // Copy over x vector
    vector<double> tmpx {x, x + cmfd::dim};

    // Loop around matrix rows
    for (int irow = 0; irow < cmfd::dim; irow++) {
      // Get index of diagonal for current row
      int didx = get_diagonal_index(irow);

      // Perform temporary sums, first do left of diag, then right of diag
      double tmp1 = 0.0;
      for (int icol = cmfd::indptr[irow]; icol < didx; icol++)
        tmp1 += A_data[icol] * x[cmfd::indices[icol]];
      for (int icol = didx + 1; icol < cmfd::indptr[irow + 1]; icol++)
        tmp1 += A_data[icol] * x[cmfd::indices[icol]];

      // Solve for new x
      double x1 = (b[irow] - tmp1) / A_data[didx];

      // Perform overrelaxation
      x[irow] = (1.0 - w) * x[irow] + w * x1;

      // Compute residual and update error
      double res = (tmpx[irow] - x[irow]) / tmpx[irow];
      err += res * res;
    }

    // Check convergence
    err = std::sqrt(err / cmfd::dim);
    if (err < tol)
      return igs;

    // Calculate new overrelaxation parameter
    w = 1.0 / (1.0 - 0.25 * cmfd::spectral * w);
  }

  // Throw error, as max iterations met
  fatal_error("Maximum Gauss-Seidel iterations encountered.");

  // Return -1 by default, although error thrown before reaching this point
  return -1;
}

//==============================================================================
// OPENMC_INITIALIZE_LINSOLVER sets the fixed variables that are used for the
// linear solver
//==============================================================================

extern "C" void openmc_initialize_linsolver(const int* indptr, int len_indptr,
  const int* indices, int n_elements, int dim, double spectral, const int* map,
  bool use_all_threads)
{
  // Store elements of indptr
  for (int i = 0; i < len_indptr; i++)
    cmfd::indptr.push_back(indptr[i]);

  // Store elements of indices
  for (int i = 0; i < n_elements; i++)
    cmfd::indices.push_back(indices[i]);

  // Set dimenion of CMFD problem and specral radius
  cmfd::dim = dim;
  cmfd::spectral = spectral;

  // Set indexmap if 1 or 2 group problem
  if (cmfd::ng == 1 || cmfd::ng == 2) {
    // Resize indexmap and set its elements
    cmfd::indexmap.resize({static_cast<size_t>(dim), 3});
    set_indexmap(map);
  }

  // Use all threads allocated to OpenMC simulation to run CMFD solver
  cmfd::use_all_threads = use_all_threads;
}

//==============================================================================
// OPENMC_RUN_LINSOLVER runs a Gauss Seidel linear solver to solve CMFD matrix
// equations
//==============================================================================

extern "C" int openmc_run_linsolver(
  const double* A_data, const double* b, double* x, double tol)
{
  switch (cmfd::ng) {
  case 1:
    return cmfd_linsolver_1g(A_data, b, x, tol);
  case 2:
    return cmfd_linsolver_2g(A_data, b, x, tol);
  default:
    return cmfd_linsolver_ng(A_data, b, x, tol);
  }
}

void free_memory_cmfd()
{
  // Clear vectors
  cmfd::indptr.clear();
  cmfd::indices.clear();
  cmfd::egrid.clear();

  // Resize tensors to be empty
  cmfd::indexmap.resize({0});

  // Set pointers to null
  cmfd::mesh = nullptr;
}

} // namespace openmc

1	#include "openmc/cmfd_solver.h"
2
3	#include <cmath>
4
5	#ifdef _OPENMP
6	#include <omp.h>
7	#endif
8	#include "openmc/tensor.h"
9
10	#include "openmc/bank.h"
11	#include "openmc/capi.h"
12	#include "openmc/constants.h"
13	#include "openmc/error.h"
14	#include "openmc/mesh.h"
15	#include "openmc/message_passing.h"
16	#include "openmc/tallies/filter_energy.h"
17	#include "openmc/tallies/filter_mesh.h"
18	#include "openmc/tallies/tally.h"
19	#include "openmc/vector.h"
20
21	namespace openmc {
22
23	namespace cmfd {
24
25	//==============================================================================
26	// Global variables
27	//==============================================================================
28
29	vector<int> indptr;
30
31	vector<int> indices;
32
33	int dim;
34
35	double spectral;
36
37	int nx, ny, nz, ng;
38
39	tensor::Tensor<int> indexmap;
40
41	int use_all_threads;
42
43	StructuredMesh* mesh;
44
45	vector<double> egrid;
46
47	double norm;
48
49	} // namespace cmfd
50
51	//==============================================================================
52	// GET_CMFD_ENERGY_BIN returns the energy bin for a source site energy
53	//==============================================================================
54
55	int get_cmfd_energy_bin(const double E)	2,250,600✔
56	{
57	// Check if energy is out of grid bounds
58	if (E < cmfd::egrid[0]) {	2,250,600!
59	// throw warning message
60	warning("Detected source point below energy grid");	×
61	return 0;	×
62	} else if (E >= cmfd::egrid[cmfd::ng]) {	2,250,600!
63	// throw warning message
64	warning("Detected source point above energy grid");	×
65	return cmfd::ng - 1;	×
66	} else {
67	// Iterate through energy grid to find matching bin
68	for (int g = 0; g < cmfd::ng; g++) {	2,571,800!
69	if (E >= cmfd::egrid[g] && E < cmfd::egrid[g + 1]) {	2,571,800!
70	return g;	2,250,600✔
71	}
72	}
73	}
74	// Return -1 by default
75	return -1;
76	}
77
78	//==============================================================================
79	// COUNT_BANK_SITES bins fission sites according to CMFD mesh and energy
80	//==============================================================================
81
82	tensor::Tensor<double> count_bank_sites(	2,299✔
83	tensor::Tensor<int>& bins, bool* outside)
84	{
85	// Determine shape of array for counts
86	std::size_t cnt_size = cmfd::nx * cmfd::ny * cmfd::nz * cmfd::ng;	2,299✔
87
88	// Create array of zeros
89	tensor::Tensor<double> cnt = tensor::zeros<double>({cnt_size});	2,299✔
90	bool outside_ = false;	2,299✔
91
92	auto bank_size = simulation::source_bank.size();	2,299✔
93	for (int i = 0; i < bank_size; i++) {	2,252,899✔
94	const auto& site = simulation::source_bank[i];	2,250,600✔
95
96	// determine scoring bin for CMFD mesh
97	int mesh_bin = cmfd::mesh->get_bin(site.r);	2,250,600✔
98
99	// if outside mesh, skip particle
100	if (mesh_bin < 0) {	2,250,600!
UNCOV 101	outside_ = true;	×
102	continue;	×
103	}
104
105	// determine scoring bin for CMFD energy
106	int energy_bin = get_cmfd_energy_bin(site.E);	2,250,600✔
107
108	// add to appropriate bin
109	cnt(mesh_bin * cmfd::ng + energy_bin) += site.wgt;	2,250,600✔
110
111	// store bin index which is used again when updating weights
112	bins[i] = mesh_bin * cmfd::ng + energy_bin;	2,250,600✔
113	}
114
115	int total = cnt.size();	2,299✔
116	tensor::Tensor<double> counts = tensor::zeros<double>({cnt_size});	2,299✔
117
118	#ifdef OPENMC_MPI
119	// collect values from all processors
120	MPI_Reduce(	836✔
121	cnt.data(), counts.data(), total, MPI_DOUBLE, MPI_SUM, 0, mpi::intracomm);	836✔
122
123	// Check if there were sites outside the mesh for any processor
124	MPI_Reduce(&outside_, outside, 1, MPI_C_BOOL, MPI_LOR, 0, mpi::intracomm);	836✔
125
126	#else
127	std::copy(cnt.data(), cnt.data() + total, counts.data());	1,463✔
128	*outside = outside_;	1,463✔
129	#endif
130
131	return counts;	2,299✔
132	}	2,299✔
133
134	//==============================================================================
135	// OPENMC_CMFD_REWEIGHT performs reweighting of particles in source bank
136	//==============================================================================
137
138	extern "C" void openmc_cmfd_reweight(	2,299✔
139	const bool feedback, const double* cmfd_src)
140	{
141	// Get size of source bank and cmfd_src
142	auto bank_size = simulation::source_bank.size();	2,299✔
143	std::size_t src_size = cmfd::nx * cmfd::ny * cmfd::nz * cmfd::ng;	2,299✔
144
145	// count bank sites for CMFD mesh, store bins in bank_bins for reweighting
146	tensor::Tensor<int> bank_bins = tensor::zeros<int>({bank_size});	2,299✔
147	bool sites_outside;	2,299✔
148	tensor::Tensor<double> sourcecounts =	2,299✔
149	count_bank_sites(bank_bins, &sites_outside);	2,299✔
150
151	// Compute CMFD weightfactors
152	tensor::Tensor<double> weightfactors = tensor::ones<double>({src_size});	2,299✔
153	if (mpi::master) {	2,299!
154	if (sites_outside) {	2,299!
UNCOV 155	fatal_error("Source sites outside of the CMFD mesh");	×
156	}
157
158	double norm = sourcecounts.sum() / cmfd::norm;	2,299✔
159	for (int i = 0; i < src_size; i++) {	25,883✔
160	if (sourcecounts[i] > 0 && cmfd_src[i] > 0) {	23,584!
161	weightfactors[i] = cmfd_src[i] * norm / sourcecounts[i];	21,912✔
162	}
163	}
164	}
165
166	if (!feedback)	2,299✔
167	return;	176✔
168
169	#ifdef OPENMC_MPI
170	// Send weightfactors to all processors
171	MPI_Bcast(weightfactors.data(), src_size, MPI_DOUBLE, 0, mpi::intracomm);	772✔
172	#endif
173
174	// Iterate through fission bank and update particle weights
175	for (int64_t i = 0; i < bank_size; i++) {	2,076,723✔
176	auto& site = simulation::source_bank[i];	2,074,600✔
177	site.wgt *= weightfactors(bank_bins(i));	2,074,600✔
178	}
179	}	6,897✔
180
181	//==============================================================================
182	// OPENMC_INITIALIZE_MESH_EGRID sets the mesh and energy grid for CMFD reweight
183	//==============================================================================
184
185	extern "C" void openmc_initialize_mesh_egrid(	165✔
186	const int meshtally_id, const int* cmfd_indices, const double norm)
187	{
188	// Make sure all CMFD memory is freed
189	free_memory_cmfd();	165✔
190
191	// Set CMFD indices
192	cmfd::nx = cmfd_indices[0];	165✔
193	cmfd::ny = cmfd_indices[1];	165✔
194	cmfd::nz = cmfd_indices[2];	165✔
195	cmfd::ng = cmfd_indices[3];	165✔
196
197	// Set CMFD reweight properties
198	cmfd::norm = norm;	165✔
199
200	// Find index corresponding to tally id
201	int32_t tally_index;	165✔
202	openmc_get_tally_index(meshtally_id, &tally_index);	165✔
203
204	// Get filters assocaited with tally
205	const auto& tally_filters = model::tallies[tally_index]->filters();	165✔
206
207	// Get mesh filter index
208	auto meshfilter_index = tally_filters[0];	165✔
209
210	// Store energy filter index if defined, otherwise set to -1
211	auto energy_index = (tally_filters.size() == 2) ? tally_filters[1] : -1;	165✔
212
213	// Get mesh index from mesh filter index
214	int32_t mesh_index;	165✔
215	openmc_mesh_filter_get_mesh(meshfilter_index, &mesh_index);	165✔
216
217	// Get mesh from mesh index
218	cmfd::mesh = dynamic_cast<StructuredMesh*>(model::meshes[mesh_index].get());	165!
219
220	// Get energy bins from energy index, otherwise use default
221	if (energy_index != -1) {	165✔
222	auto efilt_base = model::tally_filters[energy_index].get();	22!
223	auto* efilt = dynamic_cast<EnergyFilter*>(efilt_base);	22!
224	cmfd::egrid = efilt->bins();	22✔
225	} else {
226	cmfd::egrid = {0.0, INFTY};	143✔
227	}
228	}	165✔
229
230	//==============================================================================
231	// MATRIX_TO_INDICES converts a matrix index to spatial and group
232	// indices
233	//==============================================================================
234
235	void matrix_to_indices(int irow, int& g, int& i, int& j, int& k)	81,053,040✔
236	{
237	g = irow % cmfd::ng;	81,053,040✔
238	i = cmfd::indexmap(irow / cmfd::ng, 0);	81,053,040✔
239	j = cmfd::indexmap(irow / cmfd::ng, 1);	81,053,040✔
240	k = cmfd::indexmap(irow / cmfd::ng, 2);	81,053,040✔
241	}	81,053,040✔
242
243	//==============================================================================
244	// GET_DIAGONAL_INDEX returns the index in CSR index array corresponding to
245	// the diagonal element of a specified row
246	//==============================================================================
247
248	int get_diagonal_index(int row)	44,865,756✔
249	{
250	for (int j = cmfd::indptr[row]; j < cmfd::indptr[row + 1]; j++) {	90,370,920!
251	if (cmfd::indices[j] == row)	90,370,920✔
252	return j;	44,865,756✔
253	}
254
255	// Return -1 if not found
256	return -1;
257	}
258
259	//==============================================================================
260	// SET_INDEXMAP sets the elements of indexmap based on input coremap
261	//==============================================================================
262
263	void set_indexmap(const int* coremap)	154✔
264	{
265	for (int z = 0; z < cmfd::nz; z++) {	308✔
266	for (int y = 0; y < cmfd::ny; y++) {	319✔
267	for (int x = 0; x < cmfd::nx; x++) {	1,683✔
268	int idx = (z * cmfd::ny * cmfd::nx) + (y * cmfd::nx) + x;	1,518✔
269	if (coremap[idx] != CMFD_NOACCEL) {	1,518!
270	int counter = coremap[idx];	1,518✔
271	cmfd::indexmap(counter, 0) = x;	1,518✔
272	cmfd::indexmap(counter, 1) = y;	1,518✔
273	cmfd::indexmap(counter, 2) = z;	1,518✔
274	}
275	}
276	}
277	}
278	}	154✔
279
280	//==============================================================================
281	// CMFD_LINSOLVER_1G solves a one group CMFD linear system
282	//==============================================================================
283
284	int cmfd_linsolver_1g(	47,333✔
285	const double* A_data, const double* b, double* x, double tol)
286	{
287	// Set overrelaxation parameter
288	double w = 1.0;	47,333✔
289
290	// Perform Gauss-Seidel iterations
291	for (int igs = 1; igs <= 10000; igs++) {	3,699,828!
292	double err = 0.0;	3,699,828✔
293
294	// Copy over x vector
295	vector<double> tmpx {x, x + cmfd::dim};	3,699,828✔
296
297	// Perform red/black Gauss-Seidel iterations
298	for (int irb = 0; irb < 2; irb++) {	11,099,484✔
299
300	// Loop around matrix rows
301	#pragma omp parallel for reduction(+ : err) if (cmfd::use_all_threads)	4,036,506✔
302	for (int irow = 0; irow < cmfd::dim; irow++) {	38,940,250✔
303	int g, i, j, k;	35,577,100✔
304	matrix_to_indices(irow, g, i, j, k);	35,577,100✔
305
306	// Filter out black cells
307	if ((i + j + k) % 2 != irb)	35,577,100✔
308	continue;	17,788,550✔
309
310	// Get index of diagonal for current row
311	int didx = get_diagonal_index(irow);	17,788,550✔
312
313	// Perform temporary sums, first do left of diag, then right of diag
314	double tmp1 = 0.0;	17,788,550✔
315	for (int icol = cmfd::indptr[irow]; icol < didx; icol++)	33,895,525✔
316	tmp1 += A_data[icol] * x[cmfd::indices[icol]];	16,106,975✔
317	for (int icol = didx + 1; icol < cmfd::indptr[irow + 1]; icol++)	33,895,525✔
318	tmp1 += A_data[icol] * x[cmfd::indices[icol]];	16,106,975✔
319
320	// Solve for new x
321	double x1 = (b[irow] - tmp1) / A_data[didx];	17,788,550✔
322
323	// Perform overrelaxation
324	x[irow] = (1.0 - w) * x[irow] + w * x1;	17,788,550✔
325
326	// Compute residual and update error
327	double res = (tmpx[irow] - x[irow]) / tmpx[irow];	17,788,550✔
328	err += res * res;	17,788,550✔
329	}
330	}
331
332	// Check convergence
333	err = std::sqrt(err / cmfd::dim);	3,699,828✔
334	if (err < tol)	3,699,828✔
335	return igs;	47,333✔
336
337	// Calculate new overrelaxation parameter
338	w = 1.0 / (1.0 - 0.25 * cmfd::spectral * w);	3,652,495✔
339	}	3,699,828✔
340
341	// Throw error, as max iterations met
UNCOV 342	fatal_error("Maximum Gauss-Seidel iterations encountered.");	×
343
344	// Return -1 by default, although error thrown before reaching this point
345	return -1;
346	}
347
348	//==============================================================================
349	// CMFD_LINSOLVER_2G solves a two group CMFD linear system
350	//==============================================================================
351
352	int cmfd_linsolver_2g(	704✔
353	const double* A_data, const double* b, double* x, double tol)
354	{
355	// Set overrelaxation parameter
356	double w = 1.0;	704✔
357
358	// Perform Gauss-Seidel iterations
359	for (int igs = 1; igs <= 10000; igs++) {	346,980!
360	double err = 0.0;	346,980✔
361
362	// Copy over x vector
363	vector<double> tmpx {x, x + cmfd::dim};	346,980✔
364
365	// Perform red/black Gauss-Seidel iterations
366	for (int irb = 0; irb < 2; irb++) {	1,040,940✔
367
368	// Loop around matrix rows
369	#pragma omp parallel for reduction(+ : err) if (cmfd::use_all_threads)	378,490✔
370	for (int irow = 0; irow < cmfd::dim; irow += 2) {	1,577,350✔
371	int g, i, j, k;	1,261,880✔
372	matrix_to_indices(irow, g, i, j, k);	1,261,880✔
373
374	// Filter out black cells
375	if ((i + j + k) % 2 != irb)	1,261,880✔
376	continue;	630,940✔
377
378	// Get index of diagonals for current row and next row
379	int d1idx = get_diagonal_index(irow);	630,940✔
380	int d2idx = get_diagonal_index(irow + 1);	630,940✔
381
382	// Get block diagonal
383	double m11 = A_data[d1idx]; // group 1 diagonal	630,940✔
384	double m12 =	630,940✔
385	A_data[d1idx + 1]; // group 1 right of diagonal (sorted by col)	630,940✔
386	double m21 =	630,940✔
387	A_data[d2idx - 1]; // group 2 left of diagonal (sorted by col)	630,940✔
388	double m22 = A_data[d2idx]; // group 2 diagonal	630,940✔
389
390	// Analytically invert the diagonal
391	double dm = m11 * m22 - m12 * m21;	630,940✔
392	double d11 = m22 / dm;	630,940✔
393	double d12 = -m12 / dm;	630,940✔
394	double d21 = -m21 / dm;	630,940✔
395	double d22 = m11 / dm;	630,940✔
396
397	// Perform temporary sums, first do left of diag, then right of diag
398	double tmp1 = 0.0;	630,940✔
399	double tmp2 = 0.0;	630,940✔
400	for (int icol = cmfd::indptr[irow]; icol < d1idx; icol++)	1,261,880✔
401	tmp1 += A_data[icol] * x[cmfd::indices[icol]];	630,940✔
402	for (int icol = cmfd::indptr[irow + 1]; icol < d2idx - 1; icol++)	1,261,880✔
403	tmp2 += A_data[icol] * x[cmfd::indices[icol]];	630,940✔
404	for (int icol = d1idx + 2; icol < cmfd::indptr[irow + 1]; icol++)	1,261,880✔
405	tmp1 += A_data[icol] * x[cmfd::indices[icol]];	630,940✔
406	for (int icol = d2idx + 1; icol < cmfd::indptr[irow + 2]; icol++)	1,261,880✔
407	tmp2 += A_data[icol] * x[cmfd::indices[icol]];	630,940✔
408
409	// Adjust with RHS vector
410	tmp1 = b[irow] - tmp1;	630,940✔
411	tmp2 = b[irow + 1] - tmp2;	630,940✔
412
413	// Solve for new x
414	double x1 = d11 * tmp1 + d12 * tmp2;	630,940✔
415	double x2 = d21 * tmp1 + d22 * tmp2;	630,940✔
416
417	// Perform overrelaxation
418	x[irow] = (1.0 - w) * x[irow] + w * x1;	630,940✔
419	x[irow + 1] = (1.0 - w) * x[irow + 1] + w * x2;	630,940✔
420
421	// Compute residual and update error
422	double res = (tmpx[irow] - x[irow]) / tmpx[irow];	630,940✔
423	err += res * res;	630,940✔
424	}
425	}
426
427	// Check convergence
428	err = std::sqrt(err / cmfd::dim);	346,980✔
429	if (err < tol)	346,980✔
430	return igs;	704✔
431
432	// Calculate new overrelaxation parameter
433	w = 1.0 / (1.0 - 0.25 * cmfd::spectral * w);	346,276✔
434	}	346,980✔
435
436	// Throw error, as max iterations met
UNCOV 437	fatal_error("Maximum Gauss-Seidel iterations encountered.");	×
438
439	// Return -1 by default, although error thrown before reaching this point
440	return -1;
441	}
442
443	//==============================================================================
444	// CMFD_LINSOLVER_NG solves a general CMFD linear system
445	//==============================================================================
446
447	int cmfd_linsolver_ng(	484✔
448	const double* A_data, const double* b, double* x, double tol)
449	{
450	// Set overrelaxation parameter
451	double w = 1.0;	484✔
452
453	// Perform Gauss-Seidel iterations
454	for (int igs = 1; igs <= 10000; igs++) {	245,943!
455	double err = 0.0;	245,943✔
456
457	// Copy over x vector
458	vector<double> tmpx {x, x + cmfd::dim};	245,943✔
459
460	// Loop around matrix rows
461	for (int irow = 0; irow < cmfd::dim; irow++) {	3,197,259✔
462	// Get index of diagonal for current row
463	int didx = get_diagonal_index(irow);	2,951,316✔
464
465	// Perform temporary sums, first do left of diag, then right of diag
466	double tmp1 = 0.0;	2,951,316✔
467	for (int icol = cmfd::indptr[irow]; icol < didx; icol++)	8,853,948✔
468	tmp1 += A_data[icol] * x[cmfd::indices[icol]];	5,902,632✔
469	for (int icol = didx + 1; icol < cmfd::indptr[irow + 1]; icol++)	8,853,948✔
470	tmp1 += A_data[icol] * x[cmfd::indices[icol]];	5,902,632✔
471
472	// Solve for new x
473	double x1 = (b[irow] - tmp1) / A_data[didx];	2,951,316✔
474
475	// Perform overrelaxation
476	x[irow] = (1.0 - w) * x[irow] + w * x1;	2,951,316✔
477
478	// Compute residual and update error
479	double res = (tmpx[irow] - x[irow]) / tmpx[irow];	2,951,316✔
480	err += res * res;	2,951,316✔
481	}
482
483	// Check convergence
484	err = std::sqrt(err / cmfd::dim);	245,943✔
485	if (err < tol)	245,943✔
486	return igs;	484✔
487
488	// Calculate new overrelaxation parameter
489	w = 1.0 / (1.0 - 0.25 * cmfd::spectral * w);	245,459✔
490	}	245,943✔
491
492	// Throw error, as max iterations met
UNCOV 493	fatal_error("Maximum Gauss-Seidel iterations encountered.");	×
494
495	// Return -1 by default, although error thrown before reaching this point
496	return -1;
497	}
498
499	//==============================================================================
500	// OPENMC_INITIALIZE_LINSOLVER sets the fixed variables that are used for the
501	// linear solver
502	//==============================================================================
503
504	extern "C" void openmc_initialize_linsolver(const int* indptr, int len_indptr,	165✔
505	const int* indices, int n_elements, int dim, double spectral, const int* map,
506	bool use_all_threads)
507	{
508	// Store elements of indptr
509	for (int i = 0; i < len_indptr; i++)	2,024✔
510	cmfd::indptr.push_back(indptr[i]);	1,859✔
511
512	// Store elements of indices
513	for (int i = 0; i < n_elements; i++)	5,313✔
514	cmfd::indices.push_back(indices[i]);	5,148✔
515
516	// Set dimenion of CMFD problem and specral radius
517	cmfd::dim = dim;	165✔
518	cmfd::spectral = spectral;	165✔
519
520	// Set indexmap if 1 or 2 group problem
521	if (cmfd::ng == 1 \|\| cmfd::ng == 2) {	165✔
522	// Resize indexmap and set its elements
523	cmfd::indexmap.resize({static_cast<size_t>(dim), 3});	154✔
524	set_indexmap(map);	154✔
525	}
526
527	// Use all threads allocated to OpenMC simulation to run CMFD solver
528	cmfd::use_all_threads = use_all_threads;	165✔
529	}	165✔
530
531	//==============================================================================
532	// OPENMC_RUN_LINSOLVER runs a Gauss Seidel linear solver to solve CMFD matrix
533	// equations
534	//==============================================================================
535
536	extern "C" int openmc_run_linsolver(	48,521✔
537	const double* A_data, const double* b, double* x, double tol)
538	{
539	switch (cmfd::ng) {	48,521✔
540	case 1:	47,333✔
541	return cmfd_linsolver_1g(A_data, b, x, tol);	47,333✔
542	case 2:	704✔
543	return cmfd_linsolver_2g(A_data, b, x, tol);	704✔
544	default:	484✔
545	return cmfd_linsolver_ng(A_data, b, x, tol);	484✔
546	}
547	}
548
549	void free_memory_cmfd()	7,567✔
550	{
551	// Clear vectors
552	cmfd::indptr.clear();	7,567✔
553	cmfd::indices.clear();	7,567✔
554	cmfd::egrid.clear();	7,567✔
555
556	// Resize tensors to be empty
557	cmfd::indexmap.resize({0});	7,567✔
558
559	// Set pointers to null
560	cmfd::mesh = nullptr;	7,567✔
561	}	7,567✔
562
563	} // namespace openmc

openmc-dev / openmc / 23084721708

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous