From 3f0878dd367b556989ea2ab184bf2bbe519995f4 Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Thu, 1 Nov 2018 02:54:21 -0400 Subject: [PATCH] Precompute array indices, clean up readability more, fix bug with coremap axis swapping --- openmc/cmfd.py | 2878 +++++++++++++++++++++++++----------------------- 1 file changed, 1495 insertions(+), 1383 deletions(-) diff --git a/openmc/cmfd.py b/openmc/cmfd.py index 9a82d0fd83..3e97f6564d 100644 --- a/openmc/cmfd.py +++ b/openmc/cmfd.py @@ -840,6 +840,7 @@ class CMFDRun(object): self._time_cmfdbuild = None self._time_cmfdsolve = None self._intracomm = None + # TODO Add all index related variables @property def cmfd_begin(self): @@ -1052,6 +1053,12 @@ class CMFDRun(object): # Configure CMFD parameters and tallies self._configure_cmfd() + # Set up CMFD coremap + self._set_coremap() + + # TEMP + self._precompute_array_indices() + # Initialize simulation openmc.capi.simulation_init() @@ -1286,10 +1293,6 @@ class CMFDRun(object): or with traditional for loops """ - # Set up CMFD coremap - if (self._mat_dim == _CMFD_NOACCEL): - self._set_coremap() - # Calculate all cross sections based on reaction rates from last batch self._compute_xs() @@ -1474,6 +1477,7 @@ class CMFDRun(object): phi_g = phi[:,g] cmfd_flux[idx + (np.full((n,),g),)] = phi_g[self._coremap[idx]] cmfd_flux2[idx + (g,)] = phi_g[self._coremap[idx]] + print("CMFD Flux") print(np.where(cmfd_flux != cmfd_flux2)) # Compute fission source @@ -1608,6 +1612,7 @@ class CMFDRun(object): energy_bins2[idx] = 0 idx = np.where((source_energies >= energy[0]) & (source_energies <= energy[-1])) energy_bins2[idx] = ng - np.digitize(source_energies, energy) + print("Energy Bins Reweight") print(np.where(energy_bins != energy_bins2)) # Determine weight factor of each particle based on its mesh index @@ -1668,6 +1673,7 @@ class CMFDRun(object): energy_bins2[idx] = ng - 1 idx = np.where((source_energies >= energy[0]) & (source_energies <= energy[-1])) energy_bins2[idx] = np.digitize(source_energies, energy) - 1 + print("Energy bins bank sites") print(np.where(energy_bins != energy_bins2)) # Determine all unique combinations of mesh bin and energy bin, and @@ -1760,6 +1766,1433 @@ class CMFDRun(object): return loss, prod + def _precompute_array_indices(self): + nx = self._indices[0] + ny = self._indices[1] + nz = self._indices[2] + ng = self._indices[3] + + # Logical for determining whether region of interest is accelerated region + is_accel = self._coremap != _CMFD_NOACCEL + # Logical for determining whether a zero flux "albedo" b.c. should be + # applied + is_zero_flux_alb = abs(self._albedo - _ZERO_FLUX) < _TINY_BIT + x_inds, y_inds, z_inds = np.indices((nx, ny, nz)) + + # Define slice equivalent to _accel[0,:,:] + slice_x = x_inds[:1,:,:] + slice_y = y_inds[:1,:,:] + slice_z = z_inds[:1,:,:] + bndry_accel = is_accel[(slice_x, slice_y, slice_z)] + self._first_x_accel = (slice_x[bndry_accel], slice_y[bndry_accel], slice_z[bndry_accel]) + + # Define slice equivalent to is_accel[-1,:,:] + slice_x = x_inds[-1:,:,:] + slice_y = y_inds[-1:,:,:] + slice_z = z_inds[-1:,:,:] + bndry_accel = is_accel[(slice_x, slice_y, slice_z)] + self._last_x_accel = (slice_x[bndry_accel], slice_y[bndry_accel], slice_z[bndry_accel]) + + # Define slice equivalent to is_accel[:,0,:] + slice_x = x_inds[:,:1,:] + slice_y = y_inds[:,:1,:] + slice_z = z_inds[:,:1,:] + bndry_accel = is_accel[(slice_x, slice_y, slice_z)] + self._first_y_accel = (slice_x[bndry_accel], slice_y[bndry_accel], slice_z[bndry_accel]) + + # Define slice equivalent to is_accel[:,-1,:] + slice_x = x_inds[:,-1:,:] + slice_y = y_inds[:,-1:,:] + slice_z = z_inds[:,-1:,:] + bndry_accel = is_accel[(slice_x, slice_y, slice_z)] + self._last_y_accel = (slice_x[bndry_accel], slice_y[bndry_accel], slice_z[bndry_accel]) + + # Define slice equivalent to is_accel[:,:,0] + slice_x = x_inds[:,:,:1] + slice_y = y_inds[:,:,:1] + slice_z = z_inds[:,:,:1] + bndry_accel = is_accel[(x_inds, y_inds, 0)] + self._first_z_accel = (slice_x[bndry_accel], slice_y[bndry_accel], slice_z[bndry_accel]) + + # Define slice equivalent to is_accel[:,:,-1] + slice_x = x_inds[:,:,-1:] + slice_y = y_inds[:,:,-1:] + slice_z = z_inds[:,:,-1:] + bndry_accel = is_accel[(slice_x, slice_y, slice_z)] + self._last_z_accel = (slice_x[bndry_accel], slice_y[bndry_accel], slice_z[bndry_accel]) + + # Logical for whether neighboring cell to the left is reflector region + adj_reflector = np.roll(self._coremap, 1, axis=0) == _CMFD_NOACCEL + + # Define slice equivalent to is_accel[1:,:,:] & adj_reflector[1:,:,:] + slice_x = x_inds[1:,:,:] + slice_y = y_inds[1:,:,:] + slice_z = z_inds[1:,:,:] + bndry_accel = (is_accel[(slice_x, slice_y, slice_z)] & + adj_reflector[(slice_x, slice_y, slice_z)]) + self._notfirst_x_accel_adj_ref = (slice_x[bndry_accel], slice_y[bndry_accel], slice_z[bndry_accel]) + + # Define slice equivalent to is_accel[1:,:,:] & ~adj_reflector[1:,:,:] + bndry_accel = (is_accel[(slice_x, slice_y, slice_z)] & + ~adj_reflector[(slice_x, slice_y, slice_z)]) + self._notfirst_x_accel_not_adj_ref = (slice_x[bndry_accel], slice_y[bndry_accel], slice_z[bndry_accel]) + + # Logical for whether neighboring cell to the right is reflector region + adj_reflector = np.roll(self._coremap, -1, axis=0) == _CMFD_NOACCEL + + # Define slice equivalent to is_accel[:-1,:,:] & adj_reflector[:-1,:,:] + slice_x = x_inds[:-1,:,:] + slice_y = y_inds[:-1,:,:] + slice_z = z_inds[:-1,:,:] + bndry_accel = (is_accel[(slice_x, slice_y, slice_z)] & + adj_reflector[(slice_x, slice_y, slice_z)]) + self._notlast_x_accel_adj_ref = (slice_x[bndry_accel], slice_y[bndry_accel], slice_z[bndry_accel]) + + # Define slice equivalent to is_accel[:-1,:,:] & ~adj_reflector[:-1,:,:] + bndry_accel = (is_accel[(slice_x, slice_y, slice_z)] & + ~adj_reflector[(slice_x, slice_y, slice_z)]) + self._notlast_x_accel_not_adj_ref = (slice_x[bndry_accel], slice_y[bndry_accel], slice_z[bndry_accel]) + + # Logical for whether neighboring cell to the back is reflector region + adj_reflector = np.roll(self._coremap, 1, axis=1) == _CMFD_NOACCEL + + # Define slice equivalent to is_accel[:,1:,:] & adj_reflector[:,1:,:] + slice_x = x_inds[:,1:,:] + slice_y = y_inds[:,1:,:] + slice_z = z_inds[:,1:,:] + bndry_accel = (is_accel[(slice_x, slice_y, slice_z)] & + adj_reflector[(slice_x, slice_y, slice_z)]) + self._notfirst_y_accel_adj_ref = (slice_x[bndry_accel], slice_y[bndry_accel], slice_z[bndry_accel]) + + # Define slice equivalent to is_accel[:,1:,:] & ~adj_reflector[:,1:,:] + bndry_accel = (is_accel[(slice_x, slice_y, slice_z)] & + ~adj_reflector[(slice_x, slice_y, slice_z)]) + self._notfirst_y_accel_not_adj_ref = (slice_x[bndry_accel], slice_y[bndry_accel], slice_z[bndry_accel]) + + # Logical for whether neighboring cell to the front is reflector region + adj_reflector = np.roll(self._coremap, -1, axis=1) == _CMFD_NOACCEL + + # Define slice equivalent to is_accel[:,:-1,:] & adj_reflector[:,:-1,:] + slice_x = x_inds[:,:-1,:] + slice_y = y_inds[:,:-1,:] + slice_z = z_inds[:,:-1,:] + bndry_accel = (is_accel[(slice_x, slice_y, slice_z)] & + adj_reflector[(slice_x, slice_y, slice_z)]) + self._notlast_y_accel_adj_ref = (slice_x[bndry_accel], slice_y[bndry_accel], slice_z[bndry_accel]) + + # Define slice equivalent to is_accel[:,:-1,:] & ~adj_reflector[:,:-1,:] + bndry_accel = (is_accel[(slice_x, slice_y, slice_z)] & + ~adj_reflector[(slice_x, slice_y, slice_z)]) + self._notlast_y_accel_not_adj_ref = (slice_x[bndry_accel], slice_y[bndry_accel], slice_z[bndry_accel]) + + # Logical for whether neighboring cell to the bottom is reflector region + adj_reflector = np.roll(self._coremap, 1, axis=2) == _CMFD_NOACCEL + + # Define slice equivalent to is_accel[:,:,1:] & adj_reflector[:,:,1:] + slice_x = x_inds[:,:,1:] + slice_y = y_inds[:,:,1:] + slice_z = z_inds[:,:,1:] + bndry_accel = (is_accel[(slice_x, slice_y, slice_z)] & + adj_reflector[(slice_x, slice_y, slice_z)]) + self._notfirst_z_accel_adj_ref = (slice_x[bndry_accel], slice_y[bndry_accel], slice_z[bndry_accel]) + + # Define slice equivalent to is_accel[:,:,1:] & ~adj_reflector[:,:,1:] + bndry_accel = (is_accel[(slice_x, slice_y, slice_z)] & + ~adj_reflector[(slice_x, slice_y, slice_z)]) + self._notfirst_z_accel_not_adj_ref = (slice_x[bndry_accel], slice_y[bndry_accel], slice_z[bndry_accel]) + + # Logical for whether neighboring cell to the top is reflector region + adj_reflector = np.roll(self._coremap, -1, axis=2) == _CMFD_NOACCEL + + # Define slice equivalent to is_accel[:,:,:-1] & adj_reflector[:,:,:-1] + slice_x = x_inds[:,:,:-1] + slice_y = y_inds[:,:,:-1] + slice_z = z_inds[:,:,:-1] + bndry_accel = (is_accel[(slice_x, slice_y, slice_z)] & + adj_reflector[(slice_x, slice_y, slice_z)]) + self._notlast_z_accel_adj_ref = (slice_x[bndry_accel], slice_y[bndry_accel], slice_z[bndry_accel]) + + # Define slice equivalent to is_accel[:,:,:-1] & ~adj_reflector[:,:,:-1] + bndry_accel = (is_accel[(slice_x, slice_y, slice_z)] & + ~adj_reflector[(slice_x, slice_y, slice_z)]) + self._notlast_z_accel_not_adj_ref = (slice_x[bndry_accel], slice_y[bndry_accel], slice_z[bndry_accel]) + + # TEMP ALL ARRAY INDICES RELATED TO BUILDING MATRIX + # Shift coremap in all directions to determine whether leakage term + # should be defined for particular cell in matrix + coremap_shift_left = np.pad(self._coremap, ((1,0),(0,0),(0,0)), + mode='constant', constant_values=_CMFD_NOACCEL)[:-1,:,:] + + coremap_shift_right = np.pad(self._coremap, ((0,1),(0,0),(0,0)), + mode='constant', constant_values=_CMFD_NOACCEL)[1:,:,:] + + coremap_shift_back = np.pad(self._coremap, ((0,0),(1,0),(0,0)), + mode='constant', constant_values=_CMFD_NOACCEL)[:,:-1,:] + + coremap_shift_front = np.pad(self._coremap, ((0,0),(0,1),(0,0)), + mode='constant', constant_values=_CMFD_NOACCEL)[:,1:,:] + + coremap_shift_bottom = np.pad(self._coremap, ((0,0),(0,0),(1,0)), + mode='constant', constant_values=_CMFD_NOACCEL)[:,:,:-1] + + coremap_shift_top = np.pad(self._coremap, ((0,0),(0,0),(0,1)), + mode='constant', constant_values=_CMFD_NOACCEL)[:,:,1:] + + # Create empty row and column vectors to store for loss matrix + row = np.array([], dtype=int) + col = np.array([], dtype=int) + + self._is_accel = np.where(self._coremap != _CMFD_NOACCEL) + self._is_accel_neig_left = np.where((self._coremap != _CMFD_NOACCEL) & + (coremap_shift_left != _CMFD_NOACCEL)) + self._is_accel_neig_right = np.where((self._coremap != _CMFD_NOACCEL) & + (coremap_shift_right != _CMFD_NOACCEL)) + self._is_accel_neig_back = np.where((self._coremap != _CMFD_NOACCEL) & + (coremap_shift_back != _CMFD_NOACCEL)) + self._is_accel_neig_front = np.where((self._coremap != _CMFD_NOACCEL) & + (coremap_shift_front != _CMFD_NOACCEL)) + self._is_accel_neig_bottom = np.where((self._coremap != _CMFD_NOACCEL) & + (coremap_shift_bottom != _CMFD_NOACCEL)) + self._is_accel_neig_top = np.where((self._coremap != _CMFD_NOACCEL) & + (coremap_shift_top != _CMFD_NOACCEL)) + + for g in range(ng): + # Extract row and column data of regions where a cell and its + # neighbor to the left are both fuel regions + idx_x = ng * (self._coremap[self._is_accel_neig_left]) + g + idx_y = ng * (coremap_shift_left[self._is_accel_neig_left]) + g + row = np.append(row, idx_x) + col = np.append(col, idx_y) + + # Extract row and column data of regions where a cell and its + # neighbor to the right are both fuel regions + idx_x = ng * (self._coremap[self._is_accel_neig_right]) + g + idx_y = ng * (coremap_shift_right[self._is_accel_neig_right]) + g + row = np.append(row, idx_x) + col = np.append(col, idx_y) + + # Extract row and column data of regions where a cell and its + # neighbor to the back are both fuel regions + idx_x = ng * (self._coremap[self._is_accel_neig_back]) + g + idx_y = ng * (coremap_shift_back[self._is_accel_neig_back]) + g + row = np.append(row, idx_x) + col = np.append(col, idx_y) + + # Extract row and column data of regions where a cell and its + # neighbor to the front are both fuel regions + idx_x = ng * (self._coremap[self._is_accel_neig_front]) + g + idx_y = ng * (coremap_shift_front[self._is_accel_neig_front]) + g + row = np.append(row, idx_x) + col = np.append(col, idx_y) + + # Extract row and column data of regions where a cell and its + # neighbor to the bottom are both fuel regions + idx_x = ng * (self._coremap[self._is_accel_neig_bottom]) + g + idx_y = ng * (coremap_shift_bottom[self._is_accel_neig_bottom]) + g + row = np.append(row, idx_x) + col = np.append(col, idx_y) + + # Extract row and column data of regions where a cell and its + # neighbor to the top are both fuel regions + idx_x = ng * (self._coremap[self._is_accel_neig_top]) + g + idx_y = ng * (coremap_shift_top[self._is_accel_neig_top]) + g + row = np.append(row, idx_x) + col = np.append(col, idx_y) + + # Extract all regions where a cell is a fuel region + idx_x = ng * (self._coremap[self._is_accel]) + g + idx_y = idx_x + row = np.append(row, idx_x) + col = np.append(col, idx_y) + + for h in range(ng): + if h != g: + # Extract all regions where a cell is a fuel region + idx_x = ng * (self._coremap[self._is_accel]) + g + idx_y = ng * (self._coremap[self._is_accel]) + h + row = np.append(row, idx_x) + col = np.append(col, idx_y) + + # Store row and col as rows and columns of production matrix + self._loss_row = row + self._loss_col = col + + # Create empty row and column vectors to store for production matrix + row = np.array([], dtype=int) + col = np.array([], dtype=int) + + for g in range(ng): + for h in range(ng): + # Extract all regions where a cell is a fuel region + idx_x = ng * (self._coremap[self._is_accel]) + g + idx_y = ng * (self._coremap[self._is_accel]) + h + # Store rows, cols, and data to add to CSR matrix + row = np.append(row, idx_x) + col = np.append(col, idx_y) + + # Store row and col as rows and columns of production matrix + self._prod_row = row + self._prod_col = col + + + def _build_loss_matrix_vectorized_updated(self, adjoint, org_loss): + # Extract spatial and energy indices and define matrix dimension + ng = self._indices[3] + n = self._mat_dim*ng + + # Define data entries used to build csr matrix + data = np.array([]) + + dtilde_left = self._dtilde[:,:,:,:,0] + dtilde_right = self._dtilde[:,:,:,:,1] + dtilde_back = self._dtilde[:,:,:,:,2] + dtilde_front = self._dtilde[:,:,:,:,3] + dtilde_bottom = self._dtilde[:,:,:,:,4] + dtilde_top = self._dtilde[:,:,:,:,5] + dhat_left = self._dhat[:,:,:,:,0] + dhat_right = self._dhat[:,:,:,:,1] + dhat_back = self._dhat[:,:,:,:,2] + dhat_front = self._dhat[:,:,:,:,3] + dhat_bottom = self._dhat[:,:,:,:,4] + dhat_top = self._dhat[:,:,:,:,5] + + dx = self._hxyz[:,:,:,np.newaxis,0] + dy = self._hxyz[:,:,:,np.newaxis,1] + dz = self._hxyz[:,:,:,np.newaxis,2] + + # Define net leakage coefficient for each surface in each matrix element + jnet = (((dtilde_right + dhat_right) - (-1.0 * dtilde_left + dhat_left)) / dx + + ((dtilde_front + dhat_front) - (-1.0 * dtilde_back + dhat_back)) / dy + + ((dtilde_top + dhat_top) - (-1.0 * dtilde_bottom + dhat_bottom)) / dz) + + for g in range(ng): + # Define leakage terms that relate terms to their neighbors to the + # left + dtilde = self._dtilde[:,:,:,g,0][self._is_accel_neig_left] + dhat = self._dhat[:,:,:,g,0][self._is_accel_neig_left] + dx = self._hxyz[:,:,:,0][self._is_accel_neig_left] + vals = (-1.0 * dtilde - dhat) / dx + # Store data to add to CSR matrix + data = np.append(data, vals) + + # Define leakage terms that relate terms to their neighbors to the + # right + dtilde = self._dtilde[:,:,:,g,1][self._is_accel_neig_right] + dhat = self._dhat[:,:,:,g,1][self._is_accel_neig_right] + dx = self._hxyz[:,:,:,0][self._is_accel_neig_right] + vals = (-1.0 * dtilde + dhat) / dx + # Store data to add to CSR matrix + data = np.append(data, vals) + + # Define leakage terms that relate terms to their neighbors in the + # back + dtilde = self._dtilde[:,:,:,g,2][self._is_accel_neig_back] + dhat = self._dhat[:,:,:,g,2][self._is_accel_neig_back] + dy = self._hxyz[:,:,:,1][self._is_accel_neig_back] + vals = (-1.0 * dtilde - dhat) / dy + # Store data to add to CSR matrix + data = np.append(data, vals) + + # Define leakage terms that relate terms to their neighbors in the + # front + dtilde = self._dtilde[:,:,:,g,3][self._is_accel_neig_front] + dhat = self._dhat[:,:,:,g,3][self._is_accel_neig_front] + dy = self._hxyz[:,:,:,1][self._is_accel_neig_front] + vals = (-1.0 * dtilde + dhat) / dy + # Store data to add to CSR matrix + data = np.append(data, vals) + + # Define leakage terms that relate terms to their neighbors to the + # bottom + dtilde = self._dtilde[:,:,:,g,4][self._is_accel_neig_bottom] + dhat = self._dhat[:,:,:,g,4][self._is_accel_neig_bottom] + dz = self._hxyz[:,:,:,2][self._is_accel_neig_bottom] + vals = (-1.0 * dtilde - dhat) / dz + # Store data to add to CSR matrix + data = np.append(data, vals) + + # Define leakage terms that relate terms to their neighbors to the + # top + dtilde = self._dtilde[:,:,:,g,5][self._is_accel_neig_top] + dhat = self._dhat[:,:,:,g,5][self._is_accel_neig_top] + dz = self._hxyz[:,:,:,2][self._is_accel_neig_top] + vals = (-1.0 * dtilde + dhat) / dz + # Store data to add to CSR matrix + data = np.append(data, vals) + + # Define terms that relate to loss of neutrons in a cell. These + # correspond to all the diagonal entries of the loss matrix + jnet_g = jnet[:,:,:,g][self._is_accel] + total_xs = self._totalxs[:,:,:,g][self._is_accel] + scatt_xs = self._scattxs[:,:,:,g,g][self._is_accel] + vals = jnet_g + total_xs - scatt_xs + # Store data to add to CSR matrix + data = np.append(data, vals) + + # Define terms that relate to in-scattering from group to group. + # These terms relate a mesh index to all mesh indices with the same + # spatial dimensions but belong to a different energy group + for h in range(ng): + if h != g: + # Get scattering macro xs, transposed + if adjoint: + scatt_xs = self._scattxs[:,:,:,g,h][self._is_accel] + # Get scattering macro xs + else: + scatt_xs = self._scattxs[:,:,:,h,g][self._is_accel] + vals = -1.0 * scatt_xs + # Store data to add to CSR matrix + data = np.append(data, vals) + + # Create csr matrix + loss = sparse.csr_matrix((data, (self._loss_row, self._loss_col)), shape=(n, n)) + print("Loss Mat") + print(np.where(loss.toarray() != org_loss.toarray())) + return loss + + def _build_prod_matrix_vectorized_updated(self, adjoint, org_prod): + # Extract spatial and energy indices and define matrix dimension + ng = self._indices[3] + n = self._mat_dim*ng + + # Define rows, columns, and data used to build csr matrix + data = np.array([]) + + # Define terms that relate to fission production from group to group. + for g in range(ng): + for h in range(ng): + # Extract all regions where a cell is a fuel region + indices = np.where(self._coremap != _CMFD_NOACCEL) + idx_x = ng * (self._coremap[indices]) + g + idx_y = ng * (self._coremap[indices]) + h + # Get nu-fission macro xs, transposed + if adjoint: + vals = (self._nfissxs[:, :, :, g, h])[self._is_accel] + # Get nu-fission macro xs + else: + vals = (self._nfissxs[:, :, :, h, g])[self._is_accel] + # Store rows, cols, and data to add to CSR matrix + data = np.append(data, vals) + + # Create csr matrix + prod = sparse.csr_matrix((data, (self._prod_row, self._prod_col)), shape=(n, n)) + print("Prod Mat") + print(np.where(prod.toarray() != org_prod.toarray())) + return prod + + def _execute_power_iter(self, loss, prod): + """Main power iteration routine for the CMFD calculation + + Parameters + ---------- + loss : scipy.sparse.spmatrix + Sparse matrix storing elements of CMFD loss matrix + prod : scipy.sparse.spmatrix + Sparse matrix storing elements of CMFD production matrix + + Returns + ------- + phi_n : numpy.ndarray + Flux vector of CMFD problem + k_n : float + Eigenvalue of CMFD problem + dom : float + Dominance ratio of CMFD problem + + """ + # Get problem size + n = loss.shape[0] + + # Set up flux vectors, intital guess set to 1 + phi_n = np.ones((n,)) + phi_o = np.ones((n,)) + + # Set up source vectors + s_n = np.zeros((n,)) + s_o = np.zeros((n,)) + serr_v = np.zeros((n,)) + + # Set initial guess + k_n = openmc.capi.keff_temp()[0] + k_o = k_n + dw = self._cmfd_shift + k_s = k_o + dw + k_ln = 1.0/(1.0/k_n - 1.0/k_s) + k_lo = k_ln + + # Set norms to 0 + norm_n = 0.0 + norm_o = 0.0 + + # Maximum number of power iterations + maxits = 10000 + + # Perform Wielandt shift + loss -= 1.0/k_s*prod + + # Begin power iteration + for i in range(maxits): + # Check if reach max number of iterations + if i == maxits - 1: + raise OpenMCError('Reached maximum iterations in CMFD power ' + 'iteration solver.') + + # Compute source vector + s_o = prod.dot(phi_o) + + # Normalize source vector + s_o /= k_lo + + # Compute new flux vector with scipy sparse solver + phi_n = sparse.linalg.spsolve(loss, s_o) + + # Compute new source vector + s_n = prod.dot(phi_n) + + # Compute new shifted eigenvalue + k_ln = np.sum(s_n) / np.sum(s_o) + + # Compute new eigenvalue + k_n = 1.0/(1.0/k_ln + 1.0/k_s) + + # Renormalize the old source + s_o *= k_lo + + # Check convergence + iconv, norm_n = self._check_convergence(s_n, s_o, k_n, k_o, i+1) + + # If converged, calculate dominance ratio and break from loop + if iconv: + dom = norm_n / norm_o + return phi_n, k_n, dom + + # Record old values if not converged + phi_o = phi_n + k_o = k_n + k_lo = k_ln + norm_o = norm_n + + def _check_convergence(self, s_n, s_o, k_n, k_o, iter): + """Checks the convergence of the CMFD problem + + Parameters + ---------- + s_n : numpy.ndarray + Source vector from current iteration + s_o : numpy.ndarray + Source vector from previous iteration + k_n : float + K-effective from current iteration + k_o : float + K-effective from previous iteration + iter: int + Iteration number + + Returns + ------- + iconv : bool + Whether the power iteration has reached convergence + serr : float + Error in source from previous iteration to current iteration, used + for dominance ratio calculations + + """ + # Calculate error in keff + kerr = abs(k_o - k_n)/k_n + + # Calculate max error in source + serr = np.sqrt(np.sum(np.where(s_n>0, ((s_n-s_o)/s_n)**2,0))/len(s_n)) + + # Check for convergence + iconv = kerr < self._cmfd_ktol and serr < self._cmfd_stol + + # Print out to user + if self._cmfd_power_monitor and openmc.capi.master(): + str1 = ' {:d}:'.format(iter) + str2 = 'k-eff: {:0.8f}'.format(k_n) + str3 = 'k-error: {0:.5E}'.format(kerr) + str4 = 'src-error: {0:.5E}'.format(serr) + print('{0:8s}{1:20s}{2:25s}{3:s}'.format(str1, str2, str3, str4)) + sys.stdout.flush() + + return iconv, serr + + def _set_coremap(self): + """Sets the core mapping information. All regions marked with zero + are set to CMFD_NO_ACCEL, while all regions marked with 1 are set to a + unique index that maps each fuel region to a row number when building + CMFD matrices + + """ + # Set number of accelerated regions in problem. This will be related to + # the dimension of CMFD matrices + self._mat_dim = np.sum(self._coremap) + + # Define coremap as cumulative sum over accelerated regions, + # otherwise set value to _CMFD_NOACCEL + self._coremap = np.where(self._coremap==0, _CMFD_NOACCEL, + np.cumsum(self._coremap)-1) + + # Reshape coremap to three dimensional array + # Indices of coremap in user input switched in x and z axes + nx = self._indices[0] + ny = self._indices[1] + nz = self._indices[2] + self._coremap = self._coremap.reshape(nz, ny, nx) + self._coremap = np.swapaxes(self._coremap, 0, 2) + + def _compute_xs(self): + """Takes CMFD tallies from OpenMC and computes macroscopic cross + sections, flux, and diffusion coefficients for each mesh cell + + """ + # Extract energy indices + nx = self._indices[0] + ny = self._indices[1] + nz = self._indices[2] + ng = self._indices[3] + + # Set flux object and source distribution all to zeros + self._flux.fill(0.) + self._openmc_src.fill(0.) + + # Set mesh widths + self._hxyz[:,:,:,:] = openmc.capi.meshes[self._cmfd_mesh_id].width + + # Reset keff_bal to zero + self._keff_bal = 0. + + # Get tallies in-memory + tallies = openmc.capi.tallies + + # Ravel coremap as 1d array similar to how tally data is arranged + coremap = np.ravel(self._coremap.swapaxes(0, 2)) + + # Set conditional numpy array as boolean vector based on coremap + # Repeat each value for number of groups in problem + is_cmfd_accel = np.repeat(coremap != _CMFD_NOACCEL, ng) + + # Get flux from CMFD tally 0 + tally_id = self._cmfd_tally_ids[0] + tally_results = tallies[tally_id].results[:,0,1] + flux = np.where(is_cmfd_accel, tally_results, 0.) + + # Detect zero flux, abort if located + if np.any(flux[is_cmfd_accel] < _TINY_BIT): + # Get index of zero flux in flux array + idx = np.argmax(np.where(is_cmfd_accel, flux, 1) < _TINY_BIT) + + # Convert scalar idx to index in flux matrix + mat_idx = np.unravel_index(idx, self._flux.shape) + + # Throw error message (one-based indexing) + # Index of group is flipped + err_message = 'Detected zero flux without coremap overlay' + \ + ' at mesh: (' + \ + ', '.join(str(i+1) for i in mat_idx[:-1]) + \ + ') in group ' + str(ng-mat_idx[-1]) + raise OpenMCError(err_message) + + # Define target tally reshape dimensions. This defines how openmc + # tallies are ordered by dimension + target_tally_shape = [nz, ny, nx, ng] + + # Reshape flux array to target shape. Swap x and z axes so that + # flux shape is now [nx, ny, nz, ng] + reshape_flux = np.swapaxes(flux.reshape(target_tally_shape), 0, 2) + + # Flip energy axis as tally results are given in reverse order of + # energy group + self._flux = np.flip(reshape_flux, axis=3) + + # Get total rr and convert to total xs from CMFD tally 0 + tally_results = tallies[tally_id].results[:,1,1] + totalxs = np.divide(tally_results, flux, \ + where=flux>0, out=np.zeros_like(tally_results)) + + # Reshape totalxs array to target shape. Swap x and z axes so that + # shape is now [nx, ny, nz, ng] + reshape_totalxs = np.swapaxes(totalxs.reshape(target_tally_shape), 0, 2) + + # Total xs is flipped in energy axis as tally results are given in + # reverse order of energy group + self._totalxs = np.flip(reshape_totalxs, axis=3) + + # Get scattering xs from CMFD tally 1 + # flux is repeated to account for extra dimensionality of scattering xs + tally_id = self._cmfd_tally_ids[1] + tally_results = tallies[tally_id].results[:,0,1] + scattxs = np.divide(tally_results, \ + np.repeat(flux, ng), \ + where=np.repeat(flux>0, ng), \ + out=np.zeros_like(tally_results)) + + # Define target tally reshape dimensions for xs with incoming + # and outgoing energies + target_tally_shape = [nz, ny, nx, ng, ng] + + # Reshape scattxs array to target shape. Swap x and z axes so that + # shape is now [nx, ny, nz, ng, ng] + reshape_scattxs = np.swapaxes(scattxs.reshape(target_tally_shape), 0, 2) + + # Scattering xs is flipped in both incoming and outgoing energy axes + # as tally results are given in reverse order of energy group + self._scattxs = np.flip(reshape_scattxs, axis=3) + self._scattxs = np.flip(self._scattxs, axis=4) + + # Get nu-fission xs from CMFD tally 1 + # flux is repeated to account for extra dimensionality of nu-fission xs + tally_results = tallies[tally_id].results[:,1,1] + num_realizations = tallies[tally_id].num_realizations + nfissxs = np.divide(tally_results, \ + np.repeat(flux, ng), \ + where=np.repeat(flux>0, ng), \ + out=np.zeros_like(tally_results)) + + # Reshape nfissxs array to target shape. Swap x and z axes so that + # shape is now [nx, ny, nz, ng, ng] + reshape_nfissxs = np.swapaxes(nfissxs.reshape(target_tally_shape), 0, 2) + + # Nu-fission xs is flipped in both incoming and outgoing energy axes + # as tally results are given in reverse order of energy group + self._nfissxs = np.flip(reshape_nfissxs, axis=3) + self._nfissxs = np.flip(self._nfissxs, axis=4) + + # Openmc source distribution is sum of nu-fission rr in incoming energies + self._openmc_src = np.sum(self._nfissxs*self._flux[:,:,:,:,np.newaxis], + axis=3) + + # Compute k_eff from source distribution + self._keff_bal = np.sum(self._openmc_src) / num_realizations + + # Normalize openmc source distribution + self._openmc_src /= np.sum(self._openmc_src) * self._norm + + # Get surface currents from CMFD tally 2 + tally_id = self._cmfd_tally_ids[2] + tally_results = tallies[tally_id].results[:,0,1] + + # Filter tally results to include only accelerated regions + current = np.where(np.repeat(flux>0, 12), tally_results, 0.) + + # Define target tally reshape dimensions for current + target_tally_shape = [nz, ny, nx, 12, ng] + + # Reshape current array to target shape. Swap x and z axes so that + # shape is now [nx, ny, nz, ng, 12] + reshape_current = np.swapaxes(current.reshape(target_tally_shape), 0, 2) + + # Current is flipped in energy axis as tally results are given in + # reverse order of energy group + self._current = np.flip(reshape_current, axis=4) + + # Get p1 scatter xs from CMFD tally 3 + tally_id = self._cmfd_tally_ids[3] + tally_results = tallies[tally_id].results[:,0,1] + + # Define target tally reshape dimensions for p1 scatter tally + target_tally_shape = [nz, ny, nx, 2, ng] + + # Reshape and extract only p1 data from tally results (no need for p0 data) + p1scattrr = np.swapaxes(tally_results.reshape(target_tally_shape), + 0, 2)[:,:,:,1,:] + + # Store p1 scatter xs + # p1 scatter xs is flipped in energy axis as tally results are given in + # reverse order of energy group + self._p1scattxs = np.divide(np.flip(p1scattrr, axis=3), self._flux, \ + where=self._flux>0, \ + out=np.zeros_like(p1scattrr)) + + # Calculate and store diffusion coefficient + self._diffcof = np.where(self._flux>0, 1.0 / (3.0 * \ + (self._totalxs - self._p1scattxs)), 0.) + + def _compute_effective_downscatter(self): + """Changes downscatter rate for zero upscatter""" + # Extract energy index + ng = self._indices[3] + + # Return if not two groups + if ng != 2: + return + + # Extract cross sections and flux for each group + flux1 = self._flux[:,:,:,0] + flux2 = self._flux[:,:,:,1] + sigt1 = self._totalxs[:,:,:,0] + sigt2 = self._totalxs[:,:,:,1] + + # First energy index is incoming energy, second is outgoing energy + sigs11 = self._scattxs[:,:,:,0,0] + sigs21 = self._scattxs[:,:,:,1,0] + sigs12 = self._scattxs[:,:,:,0,1] + sigs22 = self._scattxs[:,:,:,1,1] + + # Compute absorption xs + siga1 = sigt1 - sigs11 - sigs12 + siga2 = sigt2 - sigs22 - sigs21 + + # Compute effective downscatter XS + sigs12_eff = sigs12 - sigs21 * np.divide(flux2, flux1, where=flux1>0, + out=np.zeros_like(flux2)) + + # Recompute total cross sections and record + self._totalxs[:,:,:,0] = siga1 + sigs11 + sigs12_eff + self._totalxs[:,:,:,1] = siga2 + sigs22 + + # Record effective dowmscatter xs + self._scattxs[:,:,:,0,1] = sigs12_eff + + # Zero out upscatter cross section + self._scattxs[:,:,:,1,0] = 0.0 + + def _neutron_balance(self): + """Computes the RMS neutron balance over the CMFD mesh""" + # Extract energy indices + ng = self._indices[3] + + # Get number of accelerated regions + num_accel = np.sum(self._coremap != _CMFD_NOACCEL) + + # Get openmc k-effective + keff = openmc.capi.keff_temp()[0] + + # Define leakage in each mesh cell and energy group + leakage = ((self._current[:,:,:,_CURRENTS['out_right'],:] - \ + self._current[:,:,:,_CURRENTS['in_right'],:]) - \ + (self._current[:,:,:,_CURRENTS['in_left'],:] - \ + self._current[:,:,:,_CURRENTS['out_left'],:])) + \ + ((self._current[:,:,:,_CURRENTS['out_front'],:] - \ + self._current[:,:,:,_CURRENTS['in_front'],:]) - \ + (self._current[:,:,:,_CURRENTS['in_back'],:] - \ + self._current[:,:,:,_CURRENTS['out_back'],:])) + \ + ((self._current[:,:,:,_CURRENTS['out_top'],:] - \ + self._current[:,:,:,_CURRENTS['in_top'],:]) - \ + (self._current[:,:,:,_CURRENTS['in_bottom'],:] - \ + self._current[:,:,:,_CURRENTS['out_bottom'],:])) + + # Compute total rr + interactions = self._totalxs * self._flux + + # Compute scattering rr by broadcasting flux in outgoing energy and + # summing over incoming energy + scattering = np.sum(self._scattxs * self._flux[:,:,:,:,np.newaxis], + axis=3) + + # Compute fission rr by broadcasting flux in outgoing energy and + # summing over incoming energy + fission = np.sum(self._nfissxs * self._flux[:,:,:,:,np.newaxis], + axis=3) + + # Compute residual + res = leakage + interactions - scattering - (1.0 / keff) * fission + + # Normalize res by flux and bank res + self._resnb = np.divide(res, self._flux, where=self._flux>0) + + # Calculate RMS and record for this batch + self._balance.append(np.sqrt( + np.sum(np.multiply(self._resnb, self._resnb)) / \ + (ng * num_accel))) + + def _compute_dtilde_vectorized_updated(self): + # Extract spatial and energy indices + nx = self._indices[0] + ny = self._indices[1] + nz = self._indices[2] + ng = self._indices[3] + self._dtilde2 = np.zeros((nx, ny, nz, ng, 6)) + + # Logical for determining whether a zero flux "albedo" b.c. should be + # applied + is_zero_flux_alb = abs(self._albedo - _ZERO_FLUX) < _TINY_BIT + + # Define dtilde at left surface for all mesh cells on left boundary + # Separate between zero flux b.c. and alebdo b.c. + boundary = self._first_x_accel + boundary_grps = boundary + (slice(None),) + D = self._diffcof[boundary_grps] + dx = self._hxyz[boundary + (np.newaxis, 0)] + if is_zero_flux_alb[0]: + self._dtilde2[boundary_grps + (0,)] = 2.0 * D / dx + else: + alb = self._albedo[0] + self._dtilde2[boundary_grps + (0,)] = ((2.0 * D * (1.0 - alb)) + / (4.0 * D * (1.0 + alb) + (1.0 - alb) * dx)) + + # Define dtilde at right surface for all mesh cells on right boundary + # Separate between zero flux b.c. and alebdo b.c. + boundary = self._last_x_accel + boundary_grps = boundary + (slice(None),) + D = self._diffcof[boundary_grps] + dx = self._hxyz[boundary + (np.newaxis, 0)] + if is_zero_flux_alb[1]: + self._dtilde2[boundary_grps + (1,)] = 2.0 * D / dx + else: + alb = self._albedo[1] + self._dtilde2[boundary_grps + (1,)] = ((2.0 * D * (1.0 - alb)) + / (4.0 * D * (1.0 + alb) + (1.0 - alb) * dx)) + + # Define dtilde at back surface for all mesh cells on back boundary + # Separate between zero flux b.c. and alebdo b.c. + boundary = self._first_y_accel + boundary_grps = boundary + (slice(None),) + D = self._diffcof[boundary_grps] + dy = self._hxyz[boundary + (np.newaxis, 1)] + if is_zero_flux_alb[2]: + self._dtilde2[boundary_grps + (2,)] = 2.0 * D / dy + else: + alb = self._albedo[2] + self._dtilde2[boundary_grps + (2,)] = ((2.0 * D * (1.0 - alb)) + / (4.0 * D * (1.0 + alb) + (1.0 - alb) * dy)) + + # Define dtilde at front surface for all mesh cells on front boundary + # Separate between zero flux b.c. and alebdo b.c. + boundary = self._last_y_accel + boundary_grps = boundary + (slice(None),) + D = self._diffcof[boundary_grps] + dy = self._hxyz[boundary + (np.newaxis, 1)] + if is_zero_flux_alb[3]: + self._dtilde2[boundary_grps + (3,)] = 2.0 * D / dy + else: + alb = self._albedo[3] + self._dtilde2[boundary_grps + (3,)] = ((2.0 * D * (1.0 - alb)) + / (4.0 * D * (1.0 + alb) + (1.0 - alb) * dy)) + + # Define dtilde at bottom surface for all mesh cells on bottom boundary + # Separate between zero flux b.c. and alebdo b.c. + boundary = self._first_z_accel + boundary_grps = boundary + (slice(None),) + D = self._diffcof[boundary_grps] + dz = self._hxyz[boundary + (np.newaxis, 2)] + if is_zero_flux_alb[4]: + self._dtilde2[boundary_grps + (4,)] = 2.0 * D / dz + else: + alb = self._albedo[4] + self._dtilde2[boundary_grps + (4,)] = ((2.0 * D * (1.0 - alb)) + / (4.0 * D * (1.0 + alb) + (1.0 - alb) * dz)) + + # Define dtilde at top surface for all mesh cells on top boundary + # Separate between zero flux b.c. and alebdo b.c. + boundary = self._last_z_accel + boundary_grps = boundary + (slice(None),) + + D = self._diffcof[boundary_grps] + dz = self._hxyz[boundary + (np.newaxis, 2)] + if is_zero_flux_alb[5]: + self._dtilde2[boundary_grps + (5,)] = 2.0 * D / dz + else: + alb = self._albedo[5] + self._dtilde2[boundary_grps + (5,)] = ((2.0 * D * (1 - alb)) + / (4.0 * D * (1.0 + alb) + (1.0 - alb) * dz)) + + # Define reflector albedo for all cells on the left surface, in case + # a cell borders a reflector region on the left + current_in_left = self._current[:,:,:,_CURRENTS['in_left'],:] + current_out_left = self._current[:,:,:,_CURRENTS['out_left'],:] + ref_albedo = np.divide(current_in_left, current_out_left, + where=current_out_left > 1.0e-10, + out=np.ones_like(current_out_left)) + + # Diffusion coefficient of neighbor to left + neig_dc = np.roll(self._diffcof, 1, axis=0) + # Cell dimensions of neighbor to left + neig_hxyz = np.roll(self._hxyz, 1, axis=0) + + # Define dtilde at left surface for all mesh cells not on left boundary + # Separate between cases where regions do and don't neighbor reflector regions + boundary = self._notfirst_x_accel_adj_ref + boundary_grps = boundary + (slice(None),) + D = self._diffcof[boundary_grps] + dx = self._hxyz[boundary + (np.newaxis, 0)] + alb = ref_albedo[boundary_grps] + self._dtilde2[boundary_grps + (0,)] = ((2.0 * D * (1.0 - alb)) + / (4.0 * D * (1.0 + alb) + (1.0 - alb) * dx)) + + boundary = self._notfirst_x_accel_not_adj_ref + boundary_grps = boundary + (slice(None),) + D = self._diffcof[boundary_grps] + dx = self._hxyz[boundary + (np.newaxis, 0)] + neig_D = neig_dc[boundary_grps] + neig_dx = neig_hxyz[boundary + (np.newaxis, 0)] + self._dtilde2[boundary_grps + (0,)] = ((2.0 * D * neig_D) + / (neig_dx * D + dx * neig_D)) + + # Define reflector albedo for all cells on the right surface, in case + # a cell borders a reflector region on the right + current_in_right = self._current[:,:,:,_CURRENTS['in_right'],:] + current_out_right = self._current[:,:,:,_CURRENTS['out_right'],:] + ref_albedo = np.divide(current_in_right, current_out_right, + where=current_out_right > 1.0e-10, + out=np.ones_like(current_out_right)) + + # Diffusion coefficient of neighbor to right + neig_dc = np.roll(self._diffcof, -1, axis=0) + # Cell dimensions of neighbor to right + neig_hxyz = np.roll(self._hxyz, -1, axis=0) + + # Define dtilde at right surface for all mesh cells not on right boundary + # Separate between cases where regions do and don't neighbor reflector regions + boundary = self._notlast_x_accel_adj_ref + boundary_grps = boundary + (slice(None),) + D = self._diffcof[boundary_grps] + dx = self._hxyz[boundary + (np.newaxis, 0)] + alb = ref_albedo[boundary_grps] + self._dtilde2[boundary_grps + (1,)] = ((2.0 * D * (1.0 - alb)) + / (4.0 * D * (1.0 + alb) + (1.0 - alb) * dx)) + + boundary = self._notlast_x_accel_not_adj_ref + boundary_grps = boundary + (slice(None),) + D = self._diffcof[boundary_grps] + dx = self._hxyz[boundary + (np.newaxis, 0)] + neig_D = neig_dc[boundary_grps] + neig_dx = neig_hxyz[boundary + (np.newaxis, 0)] + self._dtilde2[boundary_grps + (1,)] = ((2.0 * D * neig_D) + / (neig_dx * D + dx * neig_D)) + + # Define reflector albedo for all cells on the back surface, in case + # a cell borders a reflector region on the back + current_in_back = self._current[:,:,:,_CURRENTS['in_back'],:] + current_out_back = self._current[:,:,:,_CURRENTS['out_back'],:] + ref_albedo = np.divide(current_in_back, current_out_back, + where=current_out_back > 1.0e-10, + out=np.ones_like(current_out_back)) + + # Diffusion coefficient of neighbor to back + neig_dc = np.roll(self._diffcof, 1, axis=1) + # Cell dimensions of neighbor to back + neig_hxyz = np.roll(self._hxyz, 1, axis=1) + + # Define dtilde at back surface for all mesh cells not on back boundary + # Separate between cases where regions do and don't neighbor reflector regions + boundary = self._notfirst_y_accel_adj_ref + boundary_grps = boundary + (slice(None),) + D = self._diffcof[boundary_grps] + dy = self._hxyz[boundary + (np.newaxis, 1)] + alb = ref_albedo[boundary_grps] + self._dtilde2[boundary_grps + (2,)] = ((2.0 * D * (1.0 - alb)) + / (4.0 * D * (1.0 + alb) + (1.0 - alb) * dy)) + + boundary = self._notfirst_y_accel_not_adj_ref + boundary_grps = boundary + (slice(None),) + D = self._diffcof[boundary_grps] + dy = self._hxyz[boundary + (np.newaxis, 1)] + neig_D = neig_dc[boundary_grps] + neig_dy = neig_hxyz[boundary + (np.newaxis, 1)] + self._dtilde2[boundary_grps + (2,)] = ((2.0 * D * neig_D) + / (neig_dy * D + dy * neig_D)) + + # Define reflector albedo for all cells on the front surface, in case + # a cell borders a reflector region in the front + current_in_front = self._current[:,:,:,_CURRENTS['in_front'],:] + current_out_front = self._current[:,:,:,_CURRENTS['out_front'],:] + ref_albedo = np.divide(current_in_front, current_out_front, + where=current_out_front > 1.0e-10, + out=np.ones_like(current_out_front)) + + # Diffusion coefficient of neighbor to front + neig_dc = np.roll(self._diffcof, -1, axis=1) + # Cell dimensions of neighbor to front + neig_hxyz = np.roll(self._hxyz, -1, axis=1) + + # Define dtilde at front surface for all mesh cells not on front boundary + # Separate between cases where regions do and don't neighbor reflector regions + boundary = self._notlast_y_accel_adj_ref + boundary_grps = boundary + (slice(None),) + D = self._diffcof[boundary_grps] + dy = self._hxyz[boundary + (np.newaxis, 1)] + alb = ref_albedo[boundary_grps] + self._dtilde2[boundary_grps + (3,)] = ((2.0 * D * (1.0 - alb)) + / (4.0 * D * (1.0 + alb) + (1.0 - alb) * dy)) + + boundary = self._notlast_y_accel_not_adj_ref + boundary_grps = boundary + (slice(None),) + D = self._diffcof[boundary_grps] + dy = self._hxyz[boundary + (np.newaxis, 1)] + neig_D = neig_dc[boundary_grps] + neig_dy = neig_hxyz[boundary + (np.newaxis, 1)] + self._dtilde2[boundary_grps + (3,)] = ((2.0 * D * neig_D) + / (neig_dy * D + dy * neig_D)) + + # Define reflector albedo for all cells on the bottom surface, in case + # a cell borders a reflector region on the bottom + current_in_bottom = self._current[:,:,:,_CURRENTS['in_bottom'],:] + current_out_bottom = self._current[:,:,:,_CURRENTS['out_bottom'],:] + ref_albedo = np.divide(current_in_bottom, current_out_bottom, + where=current_out_bottom > 1.0e-10, + out=np.ones_like(current_out_bottom)) + + # Diffusion coefficient of neighbor to bottom + neig_dc = np.roll(self._diffcof, 1, axis=2) + # Cell dimensions of neighbor to bottom + neig_hxyz = np.roll(self._hxyz, 1, axis=2) + + # Define dtilde at bottom surface for all mesh cells not on bottom boundary + # Separate between cases where regions do and don't neighbor reflector regions + boundary = self._notfirst_z_accel_adj_ref + boundary_grps = boundary + (slice(None),) + D = self._diffcof[boundary_grps] + dz = self._hxyz[boundary + (np.newaxis, 2)] + alb = ref_albedo[boundary_grps] + self._dtilde2[boundary_grps + (4,)] = ((2.0 * D * (1.0 - alb)) + / (4.0 * D * (1.0 + alb) + (1.0 - alb) * dz)) + + boundary = self._notfirst_z_accel_not_adj_ref + boundary_grps = boundary + (slice(None),) + D = self._diffcof[boundary_grps] + dz = self._hxyz[boundary + (np.newaxis, 2)] + neig_D = neig_dc[boundary_grps] + neig_dz = neig_hxyz[boundary + (np.newaxis, 2)] + self._dtilde2[boundary_grps + (4,)] = ((2.0 * D * neig_D) + / (neig_dz * D + dz * neig_D)) + + # Define reflector albedo for all cells on the top surface, in case + # a cell borders a reflector region on the top + current_in_top = self._current[:,:,:,_CURRENTS['in_top'],:] + current_out_top = self._current[:,:,:,_CURRENTS['out_top'],:] + ref_albedo = np.divide(current_in_top, current_out_top, + where=current_out_top > 1.0e-10, + out=np.ones_like(current_out_top)) + + # Diffusion coefficient of neighbor to top + neig_dc = np.roll(self._diffcof, -1, axis=2) + # Cell dimensions of neighbor to top + neig_hxyz = np.roll(self._hxyz, -1, axis=2) + + # Define dtilde at top surface for all mesh cells not on top boundary + # Separate between cases where regions do and don't neighbor reflector regions + boundary = self._notlast_z_accel_adj_ref + boundary_grps = boundary + (slice(None),) + D = self._diffcof[boundary_grps] + dz = self._hxyz[boundary + (np.newaxis, 2)] + alb = ref_albedo[boundary_grps] + self._dtilde2[boundary_grps + (5,)] = ((2.0 * D * (1.0 - alb)) + / (4.0 * D * (1.0 + alb) + (1.0 - alb) * dz)) + + boundary = self._notlast_z_accel_not_adj_ref + boundary_grps = boundary + (slice(None),) + D = self._diffcof[boundary_grps] + dz = self._hxyz[boundary + (np.newaxis, 2)] + neig_D = neig_dc[boundary_grps] + neig_dz = neig_hxyz[boundary + (np.newaxis, 2)] + self._dtilde2[boundary_grps + (5,)] = ((2.0 * D * neig_D) + / (neig_dz * D + dz * neig_D)) + + def _compute_dhat_vectorized_updated(self): + nx = self._indices[0] + ny = self._indices[1] + nz = self._indices[2] + ng = self._indices[3] + self._dhat2 = np.zeros((nx, ny, nz, ng, 6)) + + # Define current in each direction + current_in_left = self._current[:,:,:,_CURRENTS['in_left'],:] + current_out_left = self._current[:,:,:,_CURRENTS['out_left'],:] + current_in_right = self._current[:,:,:,_CURRENTS['in_right'],:] + current_out_right = self._current[:,:,:,_CURRENTS['out_right'],:] + current_in_back = self._current[:,:,:,_CURRENTS['in_back'],:] + current_out_back = self._current[:,:,:,_CURRENTS['out_back'],:] + current_in_front = self._current[:,:,:,_CURRENTS['in_front'],:] + current_out_front = self._current[:,:,:,_CURRENTS['out_front'],:] + current_in_bottom = self._current[:,:,:,_CURRENTS['in_bottom'],:] + current_out_bottom = self._current[:,:,:,_CURRENTS['out_bottom'],:] + current_in_top = self._current[:,:,:,_CURRENTS['in_top'],:] + current_out_top = self._current[:,:,:,_CURRENTS['out_top'],:] + + dx = self._hxyz[:,:,:,np.newaxis,0] + dy = self._hxyz[:,:,:,np.newaxis,1] + dz = self._hxyz[:,:,:,np.newaxis,2] + dxdydz = np.prod(self._hxyz, axis=3)[:,:,:,np.newaxis] + + # Define net current on each face + net_current_left = (current_in_left - current_out_left) / dxdydz * dx + net_current_right = (current_out_right - current_in_right) / dxdydz * dx + net_current_back = (current_in_back - current_out_back) / dxdydz * dy + net_current_front = (current_out_front - current_in_front) / dxdydz * dy + net_current_bottom = (current_in_bottom - current_out_bottom) / dxdydz * dz + net_current_top = (current_out_top - current_in_top) / dxdydz * dz + + # Define flux in each cell + cell_flux = self._flux / dxdydz + # Extract indices of coremap that are accelerated + is_accel = self._coremap != _CMFD_NOACCEL + + # Define dhat at left surface for all mesh cells on left boundary + boundary = self._first_x_accel + boundary_grps = boundary + (slice(None),) + net_current = net_current_left[boundary_grps] + dtilde = self._dtilde[boundary + (slice(None), 0)] + flux = cell_flux[boundary_grps] + self._dhat2[boundary_grps + (0,)] = (net_current + dtilde * flux) / flux + + # Define dhat at right surface for all mesh cells on right boundary + boundary = self._last_x_accel + boundary_grps = boundary + (slice(None),) + net_current = net_current_right[boundary_grps] + dtilde = self._dtilde[boundary + (slice(None), 1)] + flux = cell_flux[boundary_grps] + self._dhat2[boundary_grps + (1,)] = (net_current - dtilde * flux) / flux + + # Define dhat at back surface for all mesh cells on back boundary + boundary = self._first_y_accel + boundary_grps = boundary + (slice(None),) + net_current = net_current_back[boundary_grps] + dtilde = self._dtilde[boundary + (slice(None), 2)] + flux = cell_flux[boundary_grps] + self._dhat2[boundary_grps + (2,)] = (net_current + dtilde * flux) / flux + + # Define dhat at front surface for all mesh cells on front boundary + boundary = self._last_y_accel + boundary_grps = boundary + (slice(None),) + net_current = net_current_front[boundary_grps] + dtilde = self._dtilde[boundary + (slice(None), 3)] + flux = cell_flux[boundary_grps] + self._dhat2[boundary_grps + (3,)] = (net_current - dtilde * flux) / flux + + # Define dhat at bottom surface for all mesh cells on bottom boundary + boundary = self._first_z_accel + boundary_grps = boundary + (slice(None),) + net_current = net_current_bottom[boundary_grps] + dtilde = self._dtilde[boundary + (slice(None), 4)] + flux = cell_flux[boundary_grps] + self._dhat2[boundary_grps + (4,)] = (net_current + dtilde * flux) / flux + + # Define dhat at top surface for all mesh cells on top boundary + boundary = self._last_z_accel + boundary_grps = boundary + (slice(None),) + net_current = net_current_top[boundary_grps] + dtilde = self._dtilde[boundary + (slice(None), 5)] + flux = cell_flux[boundary_grps] + self._dhat2[boundary_grps + (5,)] = (net_current - dtilde * flux) / flux + + # Cell flux of neighbor to left + neig_flux = np.roll(self._flux, 1, axis=0) / dxdydz + + # Define dhat at left surface for all mesh cells not on left boundary + # Separate between cases where regions do and don't neighbor reflector regions + boundary = self._notfirst_x_accel_adj_ref + boundary_grps = boundary + (slice(None),) + net_current = net_current_left[boundary_grps] + dtilde = self._dtilde[boundary + (slice(None), 0)] + flux = cell_flux[boundary_grps] + self._dhat2[boundary_grps + (0,)] = (net_current + dtilde * flux) / flux + + boundary = self._notfirst_x_accel_not_adj_ref + boundary_grps = boundary + (slice(None),) + net_current = net_current_left[boundary_grps] + dtilde = self._dtilde[boundary + (slice(None), 0)] + flux = cell_flux[boundary_grps] + flux_left = neig_flux[boundary_grps] + self._dhat2[boundary_grps + (0,)] = ((net_current - dtilde * (flux_left - flux)) + / (flux_left + flux)) + + # Cell flux of neighbor to right + neig_flux = np.roll(self._flux, -1, axis=0) / dxdydz + + # Define dhat at right surface for all mesh cells not on right boundary + # Separate between cases where regions do and don't neighbor reflector regions + boundary = self._notlast_x_accel_adj_ref + boundary_grps = boundary + (slice(None),) + net_current = net_current_right[boundary_grps] + dtilde = self._dtilde[boundary + (slice(None), 1)] + flux = cell_flux[boundary_grps] + self._dhat2[boundary_grps + (1,)] = (net_current - dtilde * flux) / flux + + boundary = self._notlast_x_accel_not_adj_ref + boundary_grps = boundary + (slice(None),) + net_current = net_current_right[boundary_grps] + dtilde = self._dtilde[boundary + (slice(None), 1)] + flux = cell_flux[boundary_grps] + flux_right = neig_flux[boundary_grps] + self._dhat2[boundary_grps + (1,)] = ((net_current + dtilde * (flux_right - flux)) + / (flux_right + flux)) + + # Cell flux of neighbor to back + neig_flux = np.roll(self._flux, 1, axis=1) / dxdydz + + # Define dhat at back surface for all mesh cells not on back boundary + # Separate between cases where regions do and don't neighbor reflector regions + boundary = self._notfirst_y_accel_adj_ref + boundary_grps = boundary + (slice(None),) + net_current = net_current_back[boundary_grps] + dtilde = self._dtilde[boundary + (slice(None), 2)] + flux = cell_flux[boundary_grps] + self._dhat2[boundary_grps + (2,)] = (net_current + dtilde * flux) / flux + + boundary = self._notfirst_y_accel_not_adj_ref + boundary_grps = boundary + (slice(None),) + net_current = net_current_back[boundary_grps] + dtilde = self._dtilde[boundary + (slice(None), 2)] + flux = cell_flux[boundary_grps] + flux_back = neig_flux[boundary_grps] + self._dhat2[boundary_grps + (2,)] = ((net_current - dtilde * (flux_back - flux)) + / (flux_back + flux)) + + # Cell flux of neighbor to front + neig_flux = np.roll(self._flux, -1, axis=1) / dxdydz + + # Define dhat at front surface for all mesh cells not on front boundary + # Separate between cases where regions do and don't neighbor reflector regions + boundary = self._notlast_y_accel_adj_ref + boundary_grps = boundary + (slice(None),) + net_current = net_current_front[boundary_grps] + dtilde = self._dtilde[boundary + (slice(None), 3)] + flux = cell_flux[boundary_grps] + self._dhat2[boundary_grps + (3,)] = (net_current - dtilde * flux) / flux + + boundary = self._notlast_y_accel_not_adj_ref + boundary_grps = boundary + (slice(None),) + net_current = net_current_front[boundary_grps] + dtilde = self._dtilde[boundary + (slice(None), 3)] + flux = cell_flux[boundary_grps] + flux_front = neig_flux[boundary_grps] + self._dhat2[boundary_grps + (3,)] = ((net_current + dtilde * (flux_front - flux)) + / (flux_front + flux)) + + # Cell flux of neighbor to bottom + neig_flux = np.roll(self._flux, 1, axis=2) / dxdydz + + # Define dhat at bottom surface for all mesh cells not on bottom boundary + # Separate between cases where regions do and don't neighbor reflector regions + boundary = self._notfirst_z_accel_adj_ref + boundary_grps = boundary + (slice(None),) + net_current = net_current_bottom[boundary_grps] + dtilde = self._dtilde[boundary + (slice(None), 4)] + flux = cell_flux[boundary_grps] + self._dhat2[boundary_grps + (4,)] = (net_current + dtilde * flux) / flux + + boundary = self._notfirst_z_accel_not_adj_ref + boundary_grps = boundary + (slice(None),) + net_current = net_current_bottom[boundary_grps] + dtilde = self._dtilde[boundary + (slice(None), 4)] + flux = cell_flux[boundary_grps] + flux_bottom = neig_flux[boundary_grps] + self._dhat2[boundary_grps + (4,)] = ((net_current - dtilde * (flux_bottom - flux)) + / (flux_bottom + flux)) + + # Cell flux of neighbor to top + neig_flux = np.roll(self._flux, -1, axis=2) / dxdydz + + # Define dhat at top surface for all mesh cells not on top boundary + # Separate between cases where regions do and don't neighbor reflector regions + boundary = self._notlast_z_accel_adj_ref + boundary_grps = boundary + (slice(None),) + net_current = net_current_top[boundary_grps] + dtilde = self._dtilde[boundary + (slice(None), 5)] + flux = cell_flux[boundary_grps] + self._dhat2[boundary_grps + (5,)] = (net_current - dtilde * flux) / flux + + boundary = self._notlast_z_accel_not_adj_ref + boundary_grps = boundary + (slice(None),) + net_current = net_current_top[boundary_grps] + dtilde = self._dtilde[boundary + (slice(None), 5)] + flux = cell_flux[boundary_grps] + flux_top = neig_flux[boundary_grps] + self._dhat2[boundary_grps + (5,)] = ((net_current + dtilde * (flux_top - flux)) + / (flux_top + flux)) + print("Dhat") + print(np.where(self._dhat2 != self._dhat)) + + def _create_cmfd_tally(self): + """Creates all tallies in-memory that are used to solve CMFD problem""" + # Create Mesh object based on CMFDMesh, stored internally + cmfd_mesh = openmc.capi.Mesh() + # Store id of Mesh object + self._cmfd_mesh_id = cmfd_mesh.id + # Set dimension and parameters of Mesh object + cmfd_mesh.dimension = self._cmfd_mesh.dimension + cmfd_mesh.set_parameters(lower_left=self._cmfd_mesh.lower_left, + upper_right=self._cmfd_mesh.upper_right, + width=self._cmfd_mesh.width) + + # Create Mesh Filter object, stored internally + mesh_filter = openmc.capi.MeshFilter() + # Set mesh for Mesh Filter + mesh_filter.mesh = cmfd_mesh + + # Set up energy filters, if applicable + if self._energy_filters: + # Create Energy Filter object, stored internally + energy_filter = openmc.capi.EnergyFilter() + # Set bins for Energy Filter + energy_filter.bins = self._egrid + + # Create Energy Out Filter object, stored internally + energyout_filter = openmc.capi.EnergyoutFilter() + # Set bins for Energy Filter + energyout_filter.bins = self._egrid + + # Create Mesh Surface Filter object, stored internally + meshsurface_filter = openmc.capi.MeshSurfaceFilter() + # Set mesh for Mesh Surface Filter + meshsurface_filter.mesh = cmfd_mesh + + # Create Legendre Filter object, stored internally + legendre_filter = openmc.capi.LegendreFilter() + # Set order for Legendre Filter + legendre_filter.order = 1 + + # Create CMFD tallies, stored internally + n_tallies = 4 + self._cmfd_tally_ids = [] + for i in range(n_tallies): + tally = openmc.capi.Tally() + # Set nuclide bins + tally.nuclides = ['total'] + self._cmfd_tally_ids.append(tally.id) + + # Set attributes of CMFD flux, total tally + if i == 0: + # Set filters for tally + if self._energy_filters: + tally.filters = [mesh_filter, energy_filter] + else: + tally.filters = [mesh_filter] + # Set scores, type, and estimator for tally + tally.scores = ['flux', 'total'] + tally.type = 'volume' + tally.estimator = 'analog' + + # Set attributes of CMFD neutron production tally + elif i == 1: + # Set filters for tally + if self._energy_filters: + tally.filters = [mesh_filter, energy_filter, energyout_filter] + else: + tally.filters = [mesh_filter] + # Set scores, type, and estimator for tally + tally.scores = ['nu-scatter', 'nu-fission'] + tally.type = 'volume' + tally.estimator = 'analog' + + # Set attributes of CMFD surface current tally + elif i == 2: + # Set filters for tally + if self._energy_filters: + tally.filters = [meshsurface_filter, energy_filter] + else: + tally.filters = [meshsurface_filter] + # Set scores, type, and estimator for tally + tally.scores = ['current'] + tally.type = 'mesh-surface' + tally.estimator = 'analog' + + # Set attributes of CMFD P1 scatter tally + elif i == 3: + # Set filters for tally + if self._energy_filters: + tally.filters = [mesh_filter, legendre_filter, energy_filter] + else: + tally.filters = [mesh_filter, legendre_filter] + # Set scores for tally + tally.scores = ['scatter'] + tally.type = 'volume' + tally.estimator = 'analog' + + # Set all tallies to be active from beginning + tally.active = True + def _build_loss_matrix(self, adjoint): """Creates matrix representing loss of neutrons. This method uses for loops to loop over all spatial indices and energy groups for easier @@ -1871,6 +3304,64 @@ class CMFDRun(object): loss = sparse.csr_matrix(loss) return loss + def _build_prod_matrix(self, adjoint): + """Creates matrix representing production of neutrons. This method uses for + loops to loop over all spatial indices and energy groups for easier + readability. Matrix is returned in CSR format by populating all entries + with a numpy array and later converting to CSR format + + Parameters + ---------- + adjoint : bool + Whether or not to run an adjoint calculation + + Returns + ------- + loss : scipy.sparse.spmatrix + Sparse matrix storing elements of CMFD production matrix + + """ + # Extract spatial and energy indices and define matrix dimension + nx = self._indices[0] + ny = self._indices[1] + nz = self._indices[2] + ng = self._indices[3] + n = self._mat_dim*ng + + # Allocate matrix + prod = np.zeros((n, n)) + + for irow in range(n): + # Get indices for row in matrix + i,j,k,g = self._matrix_to_indices(irow, nx, ny, nz, ng) + + # Check if at a reflector + if self._coremap[i,j,k] == _CMFD_NOACCEL: + continue + + # Loop around all other groups + for h in range(ng): + # Get matrix column location + hmat_idx = self._indices_to_matrix(i,j,k, h, ng) + # Check for adjoint and bank val + if adjoint: + # Get nu-fission cross section from cell, transposed! + nfissxs = self._nfissxs[i, j, k, g, h] + else: + # Get nu-fission cross section from cell + nfissxs = self._nfissxs[i, j, k, h, g] + + # Set as value to be recorded + val = nfissxs + + # record value in matrix + prod[irow, hmat_idx] = val + + # Convert matrix to csr matrix in order to use scipy sparse solver + prod = sparse.csr_matrix(prod) + + return prod + def _build_loss_matrix_vectorized(self, adjoint): """Creates matrix representing loss of neutrons. Since the end shape of the loss matrix is known a priori, this method uses a @@ -2076,274 +3567,6 @@ class CMFDRun(object): loss = sparse.csr_matrix((data, (row, col)), shape=(n, n)) return loss - def _build_loss_matrix_vectorized_updated(self, adjoint, org_loss): - # Extract spatial and energy indices and define matrix dimension - ng = self._indices[3] - n = self._mat_dim*ng - - # Define rows, columns, and data used to build csr matrix - row = np.array([], dtype=int) - col = np.array([], dtype=int) - data = np.array([]) - - dtilde_left = self._dtilde[:,:,:,:,0] - dtilde_right = self._dtilde[:,:,:,:,1] - dtilde_back = self._dtilde[:,:,:,:,2] - dtilde_front = self._dtilde[:,:,:,:,3] - dtilde_bottom = self._dtilde[:,:,:,:,4] - dtilde_top = self._dtilde[:,:,:,:,5] - dhat_left = self._dhat[:,:,:,:,0] - dhat_right = self._dhat[:,:,:,:,1] - dhat_back = self._dhat[:,:,:,:,2] - dhat_front = self._dhat[:,:,:,:,3] - dhat_bottom = self._dhat[:,:,:,:,4] - dhat_top = self._dhat[:,:,:,:,5] - - dx = self._hxyz[:,:,:,np.newaxis,0] - dy = self._hxyz[:,:,:,np.newaxis,1] - dz = self._hxyz[:,:,:,np.newaxis,2] - - # Define net leakage coefficient for each surface in each matrix element - jnet = (((dtilde_right + dhat_right) - (-1.0 * dtilde_left + dhat_left)) / dx - + ((dtilde_front + dhat_front) - (-1.0 * dtilde_back + dhat_back)) / dy - + ((dtilde_top + dhat_top) - (-1.0 * dtilde_bottom + dhat_bottom)) / dz) - - # Shift coremap in all directions to determine whether leakage term - # should be defined for particular cell in matrix - coremap_shift_left = np.pad(self._coremap, ((1,0),(0,0),(0,0)), - mode='constant', constant_values=_CMFD_NOACCEL)[:-1,:,:] - - coremap_shift_right = np.pad(self._coremap, ((0,1),(0,0),(0,0)), - mode='constant', constant_values=_CMFD_NOACCEL)[1:,:,:] - - coremap_shift_back = np.pad(self._coremap, ((0,0),(1,0),(0,0)), - mode='constant', constant_values=_CMFD_NOACCEL)[:,:-1,:] - - coremap_shift_front = np.pad(self._coremap, ((0,0),(0,1),(0,0)), - mode='constant', constant_values=_CMFD_NOACCEL)[:,1:,:] - - coremap_shift_bottom = np.pad(self._coremap, ((0,0),(0,0),(1,0)), - mode='constant', constant_values=_CMFD_NOACCEL)[:,:,:-1] - - coremap_shift_top = np.pad(self._coremap, ((0,0),(0,0),(0,1)), - mode='constant', constant_values=_CMFD_NOACCEL)[:,:,1:] - - for g in range(ng): - # Define leakage terms that relate terms to their neighbors to the - # left - - # Extract all indices where a cell and its neighbor to the left - # are both fuel regions - indices = np.where((self._coremap != _CMFD_NOACCEL) & - (coremap_shift_left != _CMFD_NOACCEL)) - idx_x = ng * (self._coremap[indices]) + g - idx_y = ng * (coremap_shift_left[indices]) + g - # Compute leakage term associated with these regions - dtilde = self._dtilde[:,:,:,g,0][indices] - dhat = self._dhat[:,:,:,g,0][indices] - dx = self._hxyz[:,:,:,0][indices] - vals = (-1.0 * dtilde - dhat) / dx - # Store rows, cols, and data to add to CSR matrix - row = np.append(row, idx_x) - col = np.append(col, idx_y) - data = np.append(data, vals) - - # Define leakage terms that relate terms to their neighbors to the - # right - - # Extract all regions where a cell and its neighbor to the right - # are both fuel regions - indices = np.where((self._coremap != _CMFD_NOACCEL) & - (coremap_shift_right != _CMFD_NOACCEL)) - idx_x = ng * (self._coremap[indices]) + g - idx_y = ng * (coremap_shift_right[indices]) + g - # Compute leakage term associated with these regions - dtilde = self._dtilde[:,:,:,g,1][indices] - dhat = self._dhat[:,:,:,g,1][indices] - dx = self._hxyz[:,:,:,0][indices] - vals = (-1.0 * dtilde + dhat) / dx - # Store rows, cols, and data to add to CSR matrix - row = np.append(row, idx_x) - col = np.append(col, idx_y) - data = np.append(data, vals) - - # Define leakage terms that relate terms to their neighbors in the - # back - - # Extract all regions where a cell and its neighbor in the back - # are both fuel regions - indices = np.where((self._coremap != _CMFD_NOACCEL) & - (coremap_shift_back != _CMFD_NOACCEL)) - idx_x = ng * (self._coremap[indices]) + g - idx_y = ng * (coremap_shift_back[indices]) + g - # Compute leakage term associated with these regions - dtilde = self._dtilde[:,:,:,g,2][indices] - dhat = self._dhat[:,:,:,g,2][indices] - dy = self._hxyz[:,:,:,1][indices] - vals = (-1.0 * dtilde - dhat) / dy - # Store rows, cols, and data to add to CSR matrix - row = np.append(row, idx_x) - col = np.append(col, idx_y) - data = np.append(data, vals) - - # Define leakage terms that relate terms to their neighbors in the - # front - - # Extract all regions where a cell and its neighbor in the front - # are both fuel regions - indices = np.where((self._coremap != _CMFD_NOACCEL) & - (coremap_shift_front != _CMFD_NOACCEL)) - idx_x = ng * (self._coremap[indices]) + g - idx_y = ng * (coremap_shift_front[indices]) + g - # Compute leakage term associated with these regions - dtilde = self._dtilde[:,:,:,g,3][indices] - dhat = self._dhat[:,:,:,g,3][indices] - dy = self._hxyz[:,:,:,1][indices] - vals = (-1.0 * dtilde + dhat) / dy - # Store rows, cols, and data to add to CSR matrix - row = np.append(row, idx_x) - col = np.append(col, idx_y) - data = np.append(data, vals) - - # Define leakage terms that relate terms to their neighbors to the - # bottom - - # Extract all regions where a cell and its neighbor to the bottom - # are both fuel regions - indices = np.where((self._coremap != _CMFD_NOACCEL) & - (coremap_shift_bottom != _CMFD_NOACCEL)) - idx_x = ng * (self._coremap[indices]) + g - idx_y = ng * (coremap_shift_bottom[indices]) + g - # Compute leakage term associated with these regions - dtilde = self._dtilde[:,:,:,g,4][indices] - dhat = self._dhat[:,:,:,g,4][indices] - dz = self._hxyz[:,:,:,2][indices] - vals = (-1.0 * dtilde - dhat) / dz - # Store rows, cols, and data to add to CSR matrix - row = np.append(row, idx_x) - col = np.append(col, idx_y) - data = np.append(data, vals) - - # Define leakage terms that relate terms to their neighbors to the - # top - - # Extract all regions where a cell and its neighbor to the - # top are both fuel regions - indices = np.where((self._coremap != _CMFD_NOACCEL) & - (coremap_shift_top != _CMFD_NOACCEL)) - idx_x = ng * (self._coremap[indices]) + g - idx_y = ng * (coremap_shift_top[indices]) + g - # Compute leakage term associated with these regions - dtilde = self._dtilde[:,:,:,g,5][indices] - dhat = self._dhat[:,:,:,g,5][indices] - dz = self._hxyz[:,:,:,2][indices] - vals = (-1.0 * dtilde + dhat) / dz - # Store rows, cols, and data to add to CSR matrix - row = np.append(row, idx_x) - col = np.append(col, idx_y) - data = np.append(data, vals) - - # Define terms that relate to loss of neutrons in a cell. These - # correspond to all the diagonal entries of the loss matrix - - # Extract all regions where a cell is a fuel region - indices = np.where(self._coremap != _CMFD_NOACCEL) - idx_x = ng * (self._coremap[indices]) + g - idx_y = idx_x - jnet_g = jnet[:,:,:,g][indices] - total_xs = self._totalxs[:,:,:,g][indices] - scatt_xs = self._scattxs[:,:,:,g,g][indices] - vals = jnet_g + total_xs - scatt_xs - # Store rows, cols, and data to add to CSR matrix - row = np.append(row, idx_x) - col = np.append(col, idx_y) - data = np.append(data, vals) - - # Define terms that relate to in-scattering from group to group. - # These terms relate a mesh index to all mesh indices with the same - # spatial dimensions but belong to a different energy group - for h in range(ng): - if h != g: - # Extract all regions where a cell is a fuel region - indices = np.where(self._coremap != _CMFD_NOACCEL) - idx_x = ng * (self._coremap[indices]) + g - idx_y = ng * (self._coremap[indices]) + h - # Get scattering macro xs, transposed - if adjoint: - scatt_xs = self._scattxs[:,:,:,g,h][indices] - # Get scattering macro xs - else: - scatt_xs = self._scattxs[:,:,:,h,g][indices] - vals = -1.0 * scatt_xs - # Store rows, cols, and data to add to CSR matrix - row = np.append(row, idx_x) - col = np.append(col, idx_y) - data = np.append(data, vals) - - # Create csr matrix - loss = sparse.csr_matrix((data, (row, col)), shape=(n, n)) - print(np.where(loss.toarray() != org_loss.toarray())) - return loss - - - def _build_prod_matrix(self, adjoint): - """Creates matrix representing production of neutrons. This method uses for - loops to loop over all spatial indices and energy groups for easier - readability. Matrix is returned in CSR format by populating all entries - with a numpy array and later converting to CSR format - - Parameters - ---------- - adjoint : bool - Whether or not to run an adjoint calculation - - Returns - ------- - loss : scipy.sparse.spmatrix - Sparse matrix storing elements of CMFD production matrix - - """ - # Extract spatial and energy indices and define matrix dimension - nx = self._indices[0] - ny = self._indices[1] - nz = self._indices[2] - ng = self._indices[3] - n = self._mat_dim*ng - - # Allocate matrix - prod = np.zeros((n, n)) - - for irow in range(n): - # Get indices for row in matrix - i,j,k,g = self._matrix_to_indices(irow, nx, ny, nz, ng) - - # Check if at a reflector - if self._coremap[i,j,k] == _CMFD_NOACCEL: - continue - - # Loop around all other groups - for h in range(ng): - # Get matrix column location - hmat_idx = self._indices_to_matrix(i,j,k, h, ng) - # Check for adjoint and bank val - if adjoint: - # Get nu-fission cross section from cell, transposed! - nfissxs = self._nfissxs[i, j, k, g, h] - else: - # Get nu-fission cross section from cell - nfissxs = self._nfissxs[i, j, k, h, g] - - # Set as value to be recorded - val = nfissxs - - # record value in matrix - prod[irow, hmat_idx] = val - - # Convert matrix to csr matrix in order to use scipy sparse solver - prod = sparse.csr_matrix(prod) - - return prod - def _build_prod_matrix_vectorized(self, adjoint): """Creates matrix representing production of neutrons. Since the end shape of the production matrix is known a priori, this method uses a @@ -2392,39 +3615,6 @@ class CMFDRun(object): prod = sparse.csr_matrix((data, (row, col)), shape=(n, n)) return prod - def _build_prod_matrix_vectorized_updated(self, adjoint, org_prod): - # Extract spatial and energy indices and define matrix dimension - ng = self._indices[3] - n = self._mat_dim*ng - - # Define rows, columns, and data used to build csr matrix - row = np.array([], dtype=int) - col = np.array([], dtype=int) - data = np.array([]) - - # Define terms that relate to fission production from group to group. - for g in range(ng): - for h in range(ng): - # Extract all regions where a cell is a fuel region - indices = np.where(self._coremap != _CMFD_NOACCEL) - idx_x = ng * (self._coremap[indices]) + g - idx_y = ng * (self._coremap[indices]) + h - # Get nu-fission macro xs, transposed - if adjoint: - vals = (self._nfissxs[:, :, :, g, h])[indices] - # Get nu-fission macro xs - else: - vals = (self._nfissxs[:, :, :, h, g])[indices] - # Store rows, cols, and data to add to CSR matrix - row = np.append(row, idx_x) - col = np.append(col, idx_y) - data = np.append(data, vals) - - # Create csr matrix - prod = sparse.csr_matrix((data, (row, col)), shape=(n, n)) - print(np.where(prod.toarray() != org_prod.toarray())) - return prod - def _matrix_to_indices(self, irow, nx, ny, nz, ng): """Converts matrix index in CMFD matrices to spatial and group indices of actual problem, based on values from coremap @@ -2489,416 +3679,6 @@ class CMFDRun(object): matidx = ng*(self._coremap[i,j,k]) + g return matidx - def _execute_power_iter(self, loss, prod): - """Main power iteration routine for the CMFD calculation - - Parameters - ---------- - loss : scipy.sparse.spmatrix - Sparse matrix storing elements of CMFD loss matrix - prod : scipy.sparse.spmatrix - Sparse matrix storing elements of CMFD production matrix - - Returns - ------- - phi_n : numpy.ndarray - Flux vector of CMFD problem - k_n : float - Eigenvalue of CMFD problem - dom : float - Dominance ratio of CMFD problem - - """ - # Get problem size - n = loss.shape[0] - - # Set up flux vectors, intital guess set to 1 - phi_n = np.ones((n,)) - phi_o = np.ones((n,)) - - # Set up source vectors - s_n = np.zeros((n,)) - s_o = np.zeros((n,)) - serr_v = np.zeros((n,)) - - # Set initial guess - k_n = openmc.capi.keff_temp()[0] - k_o = k_n - dw = self._cmfd_shift - k_s = k_o + dw - k_ln = 1.0/(1.0/k_n - 1.0/k_s) - k_lo = k_ln - - # Set norms to 0 - norm_n = 0.0 - norm_o = 0.0 - - # Maximum number of power iterations - maxits = 10000 - - # Perform Wielandt shift - loss -= 1.0/k_s*prod - - # Begin power iteration - for i in range(maxits): - # Check if reach max number of iterations - if i == maxits - 1: - raise OpenMCError('Reached maximum iterations in CMFD power ' - 'iteration solver.') - - # Compute source vector - s_o = prod.dot(phi_o) - - # Normalize source vector - s_o /= k_lo - - # Compute new flux vector with scipy sparse solver - phi_n = sparse.linalg.spsolve(loss, s_o) - - # Compute new source vector - s_n = prod.dot(phi_n) - - # Compute new shifted eigenvalue - k_ln = np.sum(s_n) / np.sum(s_o) - - # Compute new eigenvalue - k_n = 1.0/(1.0/k_ln + 1.0/k_s) - - # Renormalize the old source - s_o *= k_lo - - # Check convergence - iconv, norm_n = self._check_convergence(s_n, s_o, k_n, k_o, i+1) - - # If converged, calculate dominance ratio and break from loop - if iconv: - dom = norm_n / norm_o - return phi_n, k_n, dom - - # Record old values if not converged - phi_o = phi_n - k_o = k_n - k_lo = k_ln - norm_o = norm_n - - def _check_convergence(self, s_n, s_o, k_n, k_o, iter): - """Checks the convergence of the CMFD problem - - Parameters - ---------- - s_n : numpy.ndarray - Source vector from current iteration - s_o : numpy.ndarray - Source vector from previous iteration - k_n : float - K-effective from current iteration - k_o : float - K-effective from previous iteration - iter: int - Iteration number - - Returns - ------- - iconv : bool - Whether the power iteration has reached convergence - serr : float - Error in source from previous iteration to current iteration, used - for dominance ratio calculations - - """ - # Calculate error in keff - kerr = abs(k_o - k_n)/k_n - - # Calculate max error in source - serr = np.sqrt(np.sum(np.where(s_n>0, ((s_n-s_o)/s_n)**2,0))/len(s_n)) - - # Check for convergence - iconv = kerr < self._cmfd_ktol and serr < self._cmfd_stol - - # Print out to user - if self._cmfd_power_monitor and openmc.capi.master(): - str1 = ' {:d}:'.format(iter) - str2 = 'k-eff: {:0.8f}'.format(k_n) - str3 = 'k-error: {0:.5E}'.format(kerr) - str4 = 'src-error: {0:.5E}'.format(serr) - print('{0:8s}{1:20s}{2:25s}{3:s}'.format(str1, str2, str3, str4)) - sys.stdout.flush() - - return iconv, serr - - def _set_coremap(self): - """Sets the core mapping information. All regions marked with zero - are set to CMFD_NO_ACCEL, while all regions marked with 1 are set to a - unique index that maps each fuel region to a row number when building - CMFD matrices - - """ - # Set number of accelerated regions in problem. This will be related to - # the dimension of CMFD matrices - self._mat_dim = np.sum(self._coremap) - - # Define coremap as cumulative sum over accelerated regions, - # otherwise set value to _CMFD_NOACCEL - self._coremap = np.where(self._coremap==0, _CMFD_NOACCEL, - np.cumsum(self._coremap)-1) - - def _compute_xs(self): - """Takes CMFD tallies from OpenMC and computes macroscopic cross - sections, flux, and diffusion coefficients for each mesh cell - - """ - # Extract energy indices - nx = self._indices[0] - ny = self._indices[1] - nz = self._indices[2] - ng = self._indices[3] - - # Set flux object and source distribution all to zeros - self._flux.fill(0.) - self._openmc_src.fill(0.) - - # Set mesh widths - self._hxyz[:,:,:,:] = openmc.capi.meshes[self._cmfd_mesh_id].width - - # Reset keff_bal to zero - self._keff_bal = 0. - - # Get tallies in-memory - tallies = openmc.capi.tallies - - # Set conditional numpy array as boolean vector based on coremap - # Repeat each value for number of groups in problem - is_cmfd_accel = np.repeat(self._coremap != _CMFD_NOACCEL, ng) - - # Get flux from CMFD tally 0 - tally_id = self._cmfd_tally_ids[0] - tally_results = tallies[tally_id].results[:,0,1] - flux = np.where(is_cmfd_accel, tally_results, 0.) - - # Detect zero flux, abort if located - if np.any(flux[is_cmfd_accel] < _TINY_BIT): - # Get index of zero flux in flux array - idx = np.argmax(np.where(is_cmfd_accel, flux, 1) < _TINY_BIT) - - # Convert scalar idx to index in flux matrix - mat_idx = np.unravel_index(idx, self._flux.shape) - - # Throw error message (one-based indexing) - # Index of group is flipped - err_message = 'Detected zero flux without coremap overlay' + \ - ' at mesh: (' + \ - ', '.join(str(i+1) for i in mat_idx[:-1]) + \ - ') in group ' + str(ng-mat_idx[-1]) - raise OpenMCError(err_message) - - # Define target tally reshape dimensions. This defines how openmc - # tallies are ordered by dimension - target_tally_shape = [nz, ny, nx, ng] - - # Reshape flux array to target shape. Swap x and z axes so that - # flux shape is now [nx, ny, nz, ng] - reshape_flux = np.swapaxes(flux.reshape(target_tally_shape), 0, 2) - - # Flip energy axis as tally results are given in reverse order of - # energy group - self._flux = np.flip(reshape_flux, axis=3) - - # Get total rr and convert to total xs from CMFD tally 0 - tally_results = tallies[tally_id].results[:,1,1] - totalxs = np.divide(tally_results, flux, \ - where=flux>0, out=np.zeros_like(tally_results)) - - # Reshape totalxs array to target shape. Swap x and z axes so that - # shape is now [nx, ny, nz, ng] - reshape_totalxs = np.swapaxes(totalxs.reshape(target_tally_shape), 0, 2) - - # Total xs is flipped in energy axis as tally results are given in - # reverse order of energy group - self._totalxs = np.flip(reshape_totalxs, axis=3) - - # Get scattering xs from CMFD tally 1 - # flux is repeated to account for extra dimensionality of scattering xs - tally_id = self._cmfd_tally_ids[1] - tally_results = tallies[tally_id].results[:,0,1] - scattxs = np.divide(tally_results, \ - np.repeat(flux, ng), \ - where=np.repeat(flux>0, ng), \ - out=np.zeros_like(tally_results)) - - # Define target tally reshape dimensions for xs with incoming - # and outgoing energies - target_tally_shape = [nz, ny, nx, ng, ng] - - # Reshape scattxs array to target shape. Swap x and z axes so that - # shape is now [nx, ny, nz, ng, ng] - reshape_scattxs = np.swapaxes(scattxs.reshape(target_tally_shape), 0, 2) - - # Scattering xs is flipped in both incoming and outgoing energy axes - # as tally results are given in reverse order of energy group - self._scattxs = np.flip(reshape_scattxs, axis=3) - self._scattxs = np.flip(self._scattxs, axis=4) - - # Get nu-fission xs from CMFD tally 1 - # flux is repeated to account for extra dimensionality of nu-fission xs - tally_results = tallies[tally_id].results[:,1,1] - num_realizations = tallies[tally_id].num_realizations - nfissxs = np.divide(tally_results, \ - np.repeat(flux, ng), \ - where=np.repeat(flux>0, ng), \ - out=np.zeros_like(tally_results)) - - # Reshape nfissxs array to target shape. Swap x and z axes so that - # shape is now [nx, ny, nz, ng, ng] - reshape_nfissxs = np.swapaxes(nfissxs.reshape(target_tally_shape), 0, 2) - - # Nu-fission xs is flipped in both incoming and outgoing energy axes - # as tally results are given in reverse order of energy group - self._nfissxs = np.flip(reshape_nfissxs, axis=3) - self._nfissxs = np.flip(self._nfissxs, axis=4) - - # Openmc source distribution is sum of nu-fission rr in incoming energies - self._openmc_src = np.sum(self._nfissxs*self._flux[:,:,:,:,np.newaxis], - axis=3) - - # Compute k_eff from source distribution - self._keff_bal = np.sum(self._openmc_src) / num_realizations - - # Normalize openmc source distribution - self._openmc_src /= np.sum(self._openmc_src) * self._norm - - # Get surface currents from CMFD tally 2 - tally_id = self._cmfd_tally_ids[2] - tally_results = tallies[tally_id].results[:,0,1] - - # Filter tally results to include only accelerated regions - current = np.where(np.repeat(flux>0, 12), tally_results, 0.) - - # Define target tally reshape dimensions for current - target_tally_shape = [nz, ny, nx, 12, ng] - - # Reshape current array to target shape. Swap x and z axes so that - # shape is now [nx, ny, nz, ng, 12] - reshape_current = np.swapaxes(current.reshape(target_tally_shape), 0, 2) - - # Current is flipped in energy axis as tally results are given in - # reverse order of energy group - self._current = np.flip(reshape_current, axis=4) - - # Get p1 scatter xs from CMFD tally 3 - tally_id = self._cmfd_tally_ids[3] - tally_results = tallies[tally_id].results[:,0,1] - - # Define target tally reshape dimensions for p1 scatter tally - target_tally_shape = [nz, ny, nx, 2, ng] - - # Reshape and extract only p1 data from tally results (no need for p0 data) - p1scattrr = np.swapaxes(tally_results.reshape(target_tally_shape), - 0, 2)[:,:,:,1,:] - - # Store p1 scatter xs - # p1 scatter xs is flipped in energy axis as tally results are given in - # reverse order of energy group - self._p1scattxs = np.divide(np.flip(p1scattrr, axis=3), self._flux, \ - where=self._flux>0, \ - out=np.zeros_like(p1scattrr)) - - # Calculate and store diffusion coefficient - self._diffcof = np.where(self._flux>0, 1.0 / (3.0 * \ - (self._totalxs - self._p1scattxs)), 0.) - - # Reshape coremap to three dimensional array as all cross section data - # has been reshaped - self._coremap = self._coremap.reshape(nz,ny,nx) - self._coremap = np.swapaxes(self._coremap, 0, 2) - - def _compute_effective_downscatter(self): - """Changes downscatter rate for zero upscatter""" - # Extract energy index - ng = self._indices[3] - - # Return if not two groups - if ng != 2: - return - - # Extract cross sections and flux for each group - flux1 = self._flux[:,:,:,0] - flux2 = self._flux[:,:,:,1] - sigt1 = self._totalxs[:,:,:,0] - sigt2 = self._totalxs[:,:,:,1] - - # First energy index is incoming energy, second is outgoing energy - sigs11 = self._scattxs[:,:,:,0,0] - sigs21 = self._scattxs[:,:,:,1,0] - sigs12 = self._scattxs[:,:,:,0,1] - sigs22 = self._scattxs[:,:,:,1,1] - - # Compute absorption xs - siga1 = sigt1 - sigs11 - sigs12 - siga2 = sigt2 - sigs22 - sigs21 - - # Compute effective downscatter XS - sigs12_eff = sigs12 - sigs21 * np.divide(flux2, flux1, where=flux1>0, - out=np.zeros_like(flux2)) - - # Recompute total cross sections and record - self._totalxs[:,:,:,0] = siga1 + sigs11 + sigs12_eff - self._totalxs[:,:,:,1] = siga2 + sigs22 - - # Record effective dowmscatter xs - self._scattxs[:,:,:,0,1] = sigs12_eff - - # Zero out upscatter cross section - self._scattxs[:,:,:,1,0] = 0.0 - - def _neutron_balance(self): - """Computes the RMS neutron balance over the CMFD mesh""" - # Extract energy indices - ng = self._indices[3] - - # Get number of accelerated regions - num_accel = np.sum(self._coremap != _CMFD_NOACCEL) - - # Get openmc k-effective - keff = openmc.capi.keff_temp()[0] - - # Define leakage in each mesh cell and energy group - leakage = ((self._current[:,:,:,_CURRENTS['out_right'],:] - \ - self._current[:,:,:,_CURRENTS['in_right'],:]) - \ - (self._current[:,:,:,_CURRENTS['in_left'],:] - \ - self._current[:,:,:,_CURRENTS['out_left'],:])) + \ - ((self._current[:,:,:,_CURRENTS['out_front'],:] - \ - self._current[:,:,:,_CURRENTS['in_front'],:]) - \ - (self._current[:,:,:,_CURRENTS['in_back'],:] - \ - self._current[:,:,:,_CURRENTS['out_back'],:])) + \ - ((self._current[:,:,:,_CURRENTS['out_top'],:] - \ - self._current[:,:,:,_CURRENTS['in_top'],:]) - \ - (self._current[:,:,:,_CURRENTS['in_bottom'],:] - \ - self._current[:,:,:,_CURRENTS['out_bottom'],:])) - - # Compute total rr - interactions = self._totalxs * self._flux - - # Compute scattering rr by broadcasting flux in outgoing energy and - # summing over incoming energy - scattering = np.sum(self._scattxs * self._flux[:,:,:,:,np.newaxis], - axis=3) - - # Compute fission rr by broadcasting flux in outgoing energy and - # summing over incoming energy - fission = np.sum(self._nfissxs * self._flux[:,:,:,:,np.newaxis], - axis=3) - - # Compute residual - res = leakage + interactions - scattering - (1.0 / keff) * fission - - # Normalize res by flux and bank res - self._resnb = np.divide(res, self._flux, where=self._flux>0) - - # Calculate RMS and record for this batch - self._balance.append(np.sqrt( - np.sum(np.multiply(self._resnb, self._resnb)) / \ - (ng * num_accel))) - def _compute_dtilde_vectorized(self): """Computes the diffusion coupling coefficient using a vectorized numpy approach. Aggregate values for the dtilde multidimensional array are @@ -3136,328 +3916,6 @@ class CMFDRun(object): (neig_hxyz[:,:,:-1,np.newaxis,2] * self._diffcof[:,:,:-1,:] + \ self._hxyz[:,:,:-1,np.newaxis,2] * neig_dc[:,:,:-1,:])), 0.0) - ################################################################################# - def _compute_dtilde_vectorized_updated(self): - nx = self._indices[0] - ny = self._indices[1] - nz = self._indices[2] - ng = self._indices[3] - self._dtilde2 = np.zeros((nx, ny, nz, ng, 6)) - - # Logical for determining whether region of interest is accelerated region - is_accel = self._coremap != _CMFD_NOACCEL - # Logical for determining whether a zero flux "albedo" b.c. should be - # applied - is_zero_flux_alb = abs(self._albedo - _ZERO_FLUX) < _TINY_BIT - - # Define dtilde at left surface for all mesh cells on left boundary - # Separate between zero flux b.c. and alebdo b.c. - mask = np.full(self._coremap.shape, False, dtype=bool) - mask[0,:,:] = True - boundary = np.where(is_accel & mask) - boundary_grps = boundary + (slice(None),) - D = self._diffcof[boundary_grps] - dx = self._hxyz[boundary + (np.newaxis, 0)] - if is_zero_flux_alb[0]: - self._dtilde2[boundary_grps + (0,)] = 2.0 * D / dx - else: - alb = self._albedo[0] - self._dtilde2[boundary_grps + (0,)] = ((2.0 * D * (1.0 - alb)) - / (4.0 * D * (1.0 + alb) + (1.0 - alb) * dx)) - - # Define dtilde at right surface for all mesh cells on right boundary - # Separate between zero flux b.c. and alebdo b.c. - mask = np.full(self._coremap.shape, False, dtype=bool) - mask[-1,:,:] = True - boundary = np.where(is_accel & mask) - boundary_grps = boundary + (slice(None),) - D = self._diffcof[boundary_grps] - dx = self._hxyz[boundary + (np.newaxis, 0)] - if is_zero_flux_alb[1]: - self._dtilde2[boundary_grps + (1,)] = 2.0 * D / dx - else: - alb = self._albedo[1] - self._dtilde2[boundary_grps + (1,)] = ((2.0 * D * (1.0 - alb)) - / (4.0 * D * (1.0 + alb) + (1.0 - alb) * dx)) - - # Define dtilde at back surface for all mesh cells on back boundary - # Separate between zero flux b.c. and alebdo b.c. - mask = np.full(self._coremap.shape, False, dtype=bool) - mask[:,0,:] = True - boundary = np.where(is_accel & mask) - boundary_grps = boundary + (slice(None),) - D = self._diffcof[boundary_grps] - dy = self._hxyz[boundary + (np.newaxis, 1)] - if is_zero_flux_alb[2]: - self._dtilde2[boundary_grps + (2,)] = 2.0 * D / dy - else: - alb = self._albedo[2] - self._dtilde2[boundary_grps + (2,)] = ((2.0 * D * (1.0 - alb)) - / (4.0 * D * (1.0 + alb) + (1.0 - alb) * dy)) - - # Define dtilde at front surface for all mesh cells on front boundary - # Separate between zero flux b.c. and alebdo b.c. - mask = np.full(self._coremap.shape, False, dtype=bool) - mask[:,-1,:] = True - boundary = np.where(is_accel & mask) - boundary_grps = boundary + (slice(None),) - D = self._diffcof[boundary_grps] - dy = self._hxyz[boundary + (np.newaxis, 1)] - if is_zero_flux_alb[3]: - self._dtilde2[boundary_grps + (3,)] = 2.0 * D / dy - else: - alb = self._albedo[3] - self._dtilde2[boundary_grps + (3,)] = ((2.0 * D * (1.0 - alb)) - / (4.0 * D * (1.0 + alb) + (1.0 - alb) * dy)) - - # Define dtilde at bottom surface for all mesh cells on bottom boundary - # Separate between zero flux b.c. and alebdo b.c. - mask = np.full(self._coremap.shape, False, dtype=bool) - mask[:,:,0] = True - boundary = np.where(is_accel & mask) - boundary_grps = boundary + (slice(None),) - D = self._diffcof[boundary_grps] - dz = self._hxyz[boundary + (np.newaxis, 2)] - if is_zero_flux_alb[4]: - self._dtilde2[boundary_grps + (4,)] = 2.0 * D / dz - else: - alb = self._albedo[4] - self._dtilde2[boundary_grps + (4,)] = ((2.0 * D * (1.0 - alb)) - / (4.0 * D * (1.0 + alb) + (1.0 - alb) * dz)) - - # Define dtilde at top surface for all mesh cells on top boundary - # Separate between zero flux b.c. and alebdo b.c. - mask = np.full(self._coremap.shape, False, dtype=bool) - mask[:,:,-1] = True - boundary = np.where(is_accel & mask) - boundary_grps = boundary + (slice(None),) - - D = self._diffcof[boundary_grps] - dz = self._hxyz[boundary + (np.newaxis, 2)] - if is_zero_flux_alb[5]: - self._dtilde2[boundary_grps + (5,)] = 2.0 * D / dz - else: - alb = self._albedo[5] - self._dtilde2[boundary_grps + (5,)] = ((2.0 * D * (1 - alb)) - / (4.0 * D * (1.0 + alb) + (1.0 - alb) * dz)) - - # Define reflector albedo for all cells on the left surface, in case - # a cell borders a reflector region on the left - current_in_left = self._current[:,:,:,_CURRENTS['in_left'],:] - current_out_left = self._current[:,:,:,_CURRENTS['out_left'],:] - ref_albedo = np.divide(current_in_left, current_out_left, - where=current_out_left > 1.0e-10, - out=np.ones_like(current_out_left)) - - # Logical for whether neighboring cell to the left is reflector region - adj_reflector = np.roll(self._coremap, 1, axis=0) == _CMFD_NOACCEL - # Diffusion coefficient of neighbor to left - neig_dc = np.roll(self._diffcof, 1, axis=0) - # Cell dimensions of neighbor to left - neig_hxyz = np.roll(self._hxyz, 1, axis=0) - - # Define dtilde at left surface for all mesh cells not on left boundary - # Separate between cases where regions do and don't neighbor reflector regions - mask = np.full(self._coremap.shape, False, dtype=bool) - mask[1:,:,:] = True - boundary = np.where(is_accel & adj_reflector & mask) - boundary_grps = boundary + (slice(None),) - D = self._diffcof[boundary_grps] - dx = self._hxyz[boundary + (np.newaxis, 0)] - alb = ref_albedo[boundary_grps] - self._dtilde2[boundary_grps + (0,)] = ((2.0 * D * (1.0 - alb)) - / (4.0 * D * (1.0 + alb) + (1.0 - alb) * dx)) - - boundary = np.where(is_accel & ~adj_reflector & mask) - boundary_grps = boundary + (slice(None),) - D = self._diffcof[boundary_grps] - dx = self._hxyz[boundary + (np.newaxis, 0)] - neig_D = neig_dc[boundary_grps] - neig_dx = neig_hxyz[boundary + (np.newaxis, 0)] - self._dtilde2[boundary_grps + (0,)] = ((2.0 * D * neig_D) - / (neig_dx * D + dx * neig_D)) - - # Define reflector albedo for all cells on the right surface, in case - # a cell borders a reflector region on the right - current_in_right = self._current[:,:,:,_CURRENTS['in_right'],:] - current_out_right = self._current[:,:,:,_CURRENTS['out_right'],:] - ref_albedo = np.divide(current_in_right, current_out_right, - where=current_out_right > 1.0e-10, - out=np.ones_like(current_out_right)) - - # Logical for whether neighboring cell to the right is reflector region - adj_reflector = np.roll(self._coremap, -1, axis=0) == _CMFD_NOACCEL - # Diffusion coefficient of neighbor to right - neig_dc = np.roll(self._diffcof, -1, axis=0) - # Cell dimensions of neighbor to right - neig_hxyz = np.roll(self._hxyz, -1, axis=0) - - # Define dtilde at right surface for all mesh cells not on right boundary - # Separate between cases where regions do and don't neighbor reflector regions - mask = np.full(self._coremap.shape, False, dtype=bool) - mask[:-1,:,:] = True - boundary = np.where(is_accel & adj_reflector & mask) - boundary_grps = boundary + (slice(None),) - D = self._diffcof[boundary_grps] - dx = self._hxyz[boundary + (np.newaxis, 0)] - alb = ref_albedo[boundary_grps] - self._dtilde2[boundary_grps + (1,)] = ((2.0 * D * (1.0 - alb)) - / (4.0 * D * (1.0 + alb) + (1.0 - alb) * dx)) - - boundary = np.where(is_accel & ~adj_reflector & mask) - boundary_grps = boundary + (slice(None),) - D = self._diffcof[boundary_grps] - dx = self._hxyz[boundary + (np.newaxis, 0)] - neig_D = neig_dc[boundary_grps] - neig_dx = neig_hxyz[boundary + (np.newaxis, 0)] - self._dtilde2[boundary_grps + (1,)] = ((2.0 * D * neig_D) - / (neig_dx * D + dx * neig_D)) - - # Define reflector albedo for all cells on the back surface, in case - # a cell borders a reflector region on the back - current_in_back = self._current[:,:,:,_CURRENTS['in_back'],:] - current_out_back = self._current[:,:,:,_CURRENTS['out_back'],:] - ref_albedo = np.divide(current_in_back, current_out_back, - where=current_out_back > 1.0e-10, - out=np.ones_like(current_out_back)) - - # Logical for whether neighboring cell to the back is reflector region - adj_reflector = np.roll(self._coremap, 1, axis=1) == _CMFD_NOACCEL - # Diffusion coefficient of neighbor to back - neig_dc = np.roll(self._diffcof, 1, axis=1) - # Cell dimensions of neighbor to back - neig_hxyz = np.roll(self._hxyz, 1, axis=1) - - # Define dtilde at back surface for all mesh cells not on back boundary - # Separate between cases where regions do and don't neighbor reflector regions - mask = np.full(self._coremap.shape, False, dtype=bool) - mask[:,1:,:] = True - boundary = np.where(is_accel & adj_reflector & mask) - boundary_grps = boundary + (slice(None),) - D = self._diffcof[boundary_grps] - dy = self._hxyz[boundary + (np.newaxis, 1)] - alb = ref_albedo[boundary_grps] - self._dtilde2[boundary_grps + (2,)] = ((2.0 * D * (1.0 - alb)) - / (4.0 * D * (1.0 + alb) + (1.0 - alb) * dy)) - - boundary = np.where(is_accel & ~adj_reflector & mask) - boundary_grps = boundary + (slice(None),) - D = self._diffcof[boundary_grps] - dy = self._hxyz[boundary + (np.newaxis, 1)] - neig_D = neig_dc[boundary_grps] - neig_dy = neig_hxyz[boundary + (np.newaxis, 1)] - self._dtilde2[boundary_grps + (2,)] = ((2.0 * D * neig_D) - / (neig_dy * D + dy * neig_D)) - - # Define reflector albedo for all cells on the front surface, in case - # a cell borders a reflector region in the front - current_in_front = self._current[:,:,:,_CURRENTS['in_front'],:] - current_out_front = self._current[:,:,:,_CURRENTS['out_front'],:] - ref_albedo = np.divide(current_in_front, current_out_front, - where=current_out_front > 1.0e-10, - out=np.ones_like(current_out_front)) - - # Logical for whether neighboring cell to the front is reflector region - adj_reflector = np.roll(self._coremap, -1, axis=1) == _CMFD_NOACCEL - # Diffusion coefficient of neighbor to front - neig_dc = np.roll(self._diffcof, -1, axis=1) - # Cell dimensions of neighbor to front - neig_hxyz = np.roll(self._hxyz, -1, axis=1) - - # Define dtilde at front surface for all mesh cells not on front boundary - # Separate between cases where regions do and don't neighbor reflector regions - mask = np.full(self._coremap.shape, False, dtype=bool) - mask[:,:-1,:] = True - boundary = np.where(is_accel & adj_reflector & mask) - boundary_grps = boundary + (slice(None),) - D = self._diffcof[boundary_grps] - dy = self._hxyz[boundary + (np.newaxis, 1)] - alb = ref_albedo[boundary_grps] - self._dtilde2[boundary_grps + (3,)] = ((2.0 * D * (1.0 - alb)) - / (4.0 * D * (1.0 + alb) + (1.0 - alb) * dy)) - - boundary = np.where(is_accel & ~adj_reflector & mask) - boundary_grps = boundary + (slice(None),) - D = self._diffcof[boundary_grps] - dy = self._hxyz[boundary + (np.newaxis, 1)] - neig_D = neig_dc[boundary_grps] - neig_dy = neig_hxyz[boundary + (np.newaxis, 1)] - self._dtilde2[boundary_grps + (3,)] = ((2.0 * D * neig_D) - / (neig_dy * D + dy * neig_D)) - - # Define reflector albedo for all cells on the bottom surface, in case - # a cell borders a reflector region on the bottom - current_in_bottom = self._current[:,:,:,_CURRENTS['in_bottom'],:] - current_out_bottom = self._current[:,:,:,_CURRENTS['out_bottom'],:] - ref_albedo = np.divide(current_in_bottom, current_out_bottom, - where=current_out_bottom > 1.0e-10, - out=np.ones_like(current_out_bottom)) - - # Logical for whether neighboring cell to the bottom is reflector region - adj_reflector = np.roll(self._coremap, 1, axis=2) == _CMFD_NOACCEL - # Diffusion coefficient of neighbor to bottom - neig_dc = np.roll(self._diffcof, 1, axis=2) - # Cell dimensions of neighbor to bottom - neig_hxyz = np.roll(self._hxyz, 1, axis=2) - - # Define dtilde at bottom surface for all mesh cells not on bottom boundary - # Separate between cases where regions do and don't neighbor reflector regions - mask = np.full(self._coremap.shape, False, dtype=bool) - mask[:,:,1:] = True - boundary = np.where(is_accel & adj_reflector & mask) - boundary_grps = boundary + (slice(None),) - D = self._diffcof[boundary_grps] - dz = self._hxyz[boundary + (np.newaxis, 2)] - alb = ref_albedo[boundary_grps] - self._dtilde2[boundary_grps + (4,)] = ((2.0 * D * (1.0 - alb)) - / (4.0 * D * (1.0 + alb) + (1.0 - alb) * dz)) - - boundary = np.where(is_accel & ~adj_reflector & mask) - boundary_grps = boundary + (slice(None),) - D = self._diffcof[boundary_grps] - dz = self._hxyz[boundary + (np.newaxis, 2)] - neig_D = neig_dc[boundary_grps] - neig_dz = neig_hxyz[boundary + (np.newaxis, 2)] - self._dtilde2[boundary_grps + (4,)] = ((2.0 * D * neig_D) - / (neig_dz * D + dz * neig_D)) - - # Define reflector albedo for all cells on the top surface, in case - # a cell borders a reflector region on the top - current_in_top = self._current[:,:,:,_CURRENTS['in_top'],:] - current_out_top = self._current[:,:,:,_CURRENTS['out_top'],:] - ref_albedo = np.divide(current_in_top, current_out_top, - where=current_out_top > 1.0e-10, - out=np.ones_like(current_out_top)) - - # Logical for whether neighboring cell to the top is reflector region - adj_reflector = np.roll(self._coremap, -1, axis=2) == _CMFD_NOACCEL - # Diffusion coefficient of neighbor to top - neig_dc = np.roll(self._diffcof, -1, axis=2) - # Cell dimensions of neighbor to top - neig_hxyz = np.roll(self._hxyz, -1, axis=2) - - # Define dtilde at top surface for all mesh cells not on top boundary - # Separate between cases where regions do and don't neighbor reflector regions - mask = np.full(self._coremap.shape, False, dtype=bool) - mask[:,:,:-1] = True - boundary = np.where(is_accel & adj_reflector & mask) - boundary_grps = boundary + (slice(None),) - D = self._diffcof[boundary_grps] - dz = self._hxyz[boundary + (np.newaxis, 2)] - alb = ref_albedo[boundary_grps] - self._dtilde2[boundary_grps + (5,)] = ((2.0 * D * (1.0 - alb)) - / (4.0 * D * (1.0 + alb) + (1.0 - alb) * dz)) - - boundary = np.where(is_accel & ~adj_reflector & mask) - boundary_grps = boundary + (slice(None),) - D = self._diffcof[boundary_grps] - dz = self._hxyz[boundary + (np.newaxis, 2)] - neig_D = neig_dc[boundary_grps] - neig_dz = neig_hxyz[boundary + (np.newaxis, 2)] - self._dtilde2[boundary_grps + (5,)] = ((2.0 * D * neig_D) - / (neig_dz * D + dz * neig_D)) - print(np.where(self._dtilde2 != self._dtilde)) - def _compute_dtilde(self): """Computes the diffusion coupling coefficient by looping over all spatial regions and energy groups @@ -3683,255 +4141,6 @@ class CMFDRun(object): (neig_flux[:,:,:-1,:] - cell_flux[:,:,:-1,:])) / \ (neig_flux[:,:,:-1,:] + cell_flux[:,:,:-1,:])), 0.0) - ################################################################################# - def _compute_dhat_vectorized_updated(self): - nx = self._indices[0] - ny = self._indices[1] - nz = self._indices[2] - ng = self._indices[3] - self._dhat2 = np.zeros((nx, ny, nz, ng, 6)) - - # Define current in each direction - current_in_left = self._current[:,:,:,_CURRENTS['in_left'],:] - current_out_left = self._current[:,:,:,_CURRENTS['out_left'],:] - current_in_right = self._current[:,:,:,_CURRENTS['in_right'],:] - current_out_right = self._current[:,:,:,_CURRENTS['out_right'],:] - current_in_back = self._current[:,:,:,_CURRENTS['in_back'],:] - current_out_back = self._current[:,:,:,_CURRENTS['out_back'],:] - current_in_front = self._current[:,:,:,_CURRENTS['in_front'],:] - current_out_front = self._current[:,:,:,_CURRENTS['out_front'],:] - current_in_bottom = self._current[:,:,:,_CURRENTS['in_bottom'],:] - current_out_bottom = self._current[:,:,:,_CURRENTS['out_bottom'],:] - current_in_top = self._current[:,:,:,_CURRENTS['in_top'],:] - current_out_top = self._current[:,:,:,_CURRENTS['out_top'],:] - - dx = self._hxyz[:,:,:,np.newaxis,0] - dy = self._hxyz[:,:,:,np.newaxis,1] - dz = self._hxyz[:,:,:,np.newaxis,2] - dxdydz = np.prod(self._hxyz, axis=3)[:,:,:,np.newaxis] - - # Define net current on each face - net_current_left = (current_in_left - current_out_left) / dxdydz * dx - net_current_right = (current_out_right - current_in_right) / dxdydz * dx - net_current_back = (current_in_back - current_out_back) / dxdydz * dy - net_current_front = (current_out_front - current_in_front) / dxdydz * dy - net_current_bottom = (current_in_bottom - current_out_bottom) / dxdydz * dz - net_current_top = (current_out_top - current_in_top) / dxdydz * dz - - # Define flux in each cell - cell_flux = self._flux / dxdydz - # Extract indices of coremap that are accelerated - is_accel = self._coremap != _CMFD_NOACCEL - - # Define dhat at left surface for all mesh cells on left boundary - mask = np.full(self._coremap.shape, False, dtype=bool) - mask[0,:,:] = True - boundary = np.where(is_accel & mask) - boundary_grps = boundary + (slice(None),) - net_current = net_current_left[boundary_grps] - dtilde = self._dtilde[boundary + (slice(None), 0)] - flux = cell_flux[boundary_grps] - self._dhat2[boundary_grps + (0,)] = (net_current + dtilde * flux) / flux - - # Define dhat at right surface for all mesh cells on right boundary - mask = np.full(self._coremap.shape, False, dtype=bool) - mask[-1,:,:] = True - boundary = np.where(is_accel & mask) - boundary_grps = boundary + (slice(None),) - net_current = net_current_right[boundary_grps] - dtilde = self._dtilde[boundary + (slice(None), 1)] - flux = cell_flux[boundary_grps] - self._dhat2[boundary_grps + (1,)] = (net_current - dtilde * flux) / flux - - # Define dhat at back surface for all mesh cells on back boundary - mask = np.full(self._coremap.shape, False, dtype=bool) - mask[:,0,:] = True - boundary = np.where(is_accel & mask) - boundary_grps = boundary + (slice(None),) - net_current = net_current_back[boundary_grps] - dtilde = self._dtilde[boundary + (slice(None), 2)] - flux = cell_flux[boundary_grps] - self._dhat2[boundary_grps + (2,)] = (net_current + dtilde * flux) / flux - - # Define dhat at front surface for all mesh cells on front boundary - mask = np.full(self._coremap.shape, False, dtype=bool) - mask[:,-1,:] = True - boundary = np.where(is_accel & mask) - boundary_grps = boundary + (slice(None),) - net_current = net_current_front[boundary_grps] - dtilde = self._dtilde[boundary + (slice(None), 3)] - flux = cell_flux[boundary_grps] - self._dhat2[boundary_grps + (3,)] = (net_current - dtilde * flux) / flux - - # Define dhat at bottom surface for all mesh cells on bottom boundary - mask = np.full(self._coremap.shape, False, dtype=bool) - mask[:,:,0] = True - boundary = np.where(is_accel & mask) - boundary_grps = boundary + (slice(None),) - net_current = net_current_bottom[boundary_grps] - dtilde = self._dtilde[boundary + (slice(None), 4)] - flux = cell_flux[boundary_grps] - self._dhat2[boundary_grps + (4,)] = (net_current + dtilde * flux) / flux - - # Define dhat at top surface for all mesh cells on top boundary - mask = np.full(self._coremap.shape, False, dtype=bool) - mask[:,:,-1] = True - boundary = np.where(is_accel & mask) - boundary_grps = boundary + (slice(None),) - net_current = net_current_top[boundary_grps] - dtilde = self._dtilde[boundary + (slice(None), 5)] - flux = cell_flux[boundary_grps] - self._dhat2[boundary_grps + (5,)] = (net_current - dtilde * flux) / flux - - # Logical for whether neighboring cell to the left is reflector region - adj_reflector = np.roll(self._coremap, 1, axis=0) == _CMFD_NOACCEL - # Cell flux of neighbor to left - neig_flux = np.roll(self._flux, 1, axis=0) / dxdydz - - # Define dhat at left surface for all mesh cells not on left boundary - # Separate between cases where regions do and don't neighbor reflector regions - mask = np.full(self._coremap.shape, False, dtype=bool) - mask[1:,:,:] = True - boundary = np.where(is_accel & adj_reflector & mask) - boundary_grps = boundary + (slice(None),) - net_current = net_current_left[boundary_grps] - dtilde = self._dtilde[boundary + (slice(None), 0)] - flux = cell_flux[boundary_grps] - self._dhat2[boundary_grps + (0,)] = (net_current + dtilde * flux) / flux - - boundary = np.where(is_accel & ~adj_reflector & mask) - boundary_grps = boundary + (slice(None),) - net_current = net_current_left[boundary_grps] - dtilde = self._dtilde[boundary + (slice(None), 0)] - flux = cell_flux[boundary_grps] - flux_left = neig_flux[boundary_grps] - self._dhat2[boundary_grps + (0,)] = ((net_current - dtilde * (flux_left - flux)) - / (flux_left + flux)) - - # Logical for whether neighboring cell to the right is reflector region - adj_reflector = np.roll(self._coremap, -1, axis=0) == _CMFD_NOACCEL - # Cell flux of neighbor to right - neig_flux = np.roll(self._flux, -1, axis=0) / dxdydz - # Define dhat at right surface for all mesh cells not on right boundary - # Separate between cases where regions do and don't neighbor reflector regions - boundary = np.where(is_accel[:-1,:,:] & adj_reflector[:-1,:,:]) - boundary_grps = boundary + (slice(None),) - net_current = net_current_right[boundary_grps] - dtilde = self._dtilde[boundary + (slice(None), 1)] - flux = cell_flux[boundary_grps] - self._dhat2[boundary_grps + (1,)] = (net_current - dtilde * flux) / flux - - boundary = np.where(is_accel[:-1,:,:] & ~adj_reflector[:-1,:,:]) - boundary_grps = boundary + (slice(None),) - net_current = net_current_right[boundary_grps] - dtilde = self._dtilde[boundary + (slice(None), 1)] - flux = cell_flux[boundary_grps] - flux_right = neig_flux[boundary_grps] - self._dhat2[boundary_grps + (1,)] = ((net_current + dtilde * (flux_right - flux)) - / (flux_right + flux)) - - # Logical for whether neighboring cell to the back is reflector region - adj_reflector = np.roll(self._coremap, 1, axis=1) == _CMFD_NOACCEL - # Cell flux of neighbor to back - neig_flux = np.roll(self._flux, 1, axis=1) / dxdydz - - # Define dhat at back surface for all mesh cells not on back boundary - # Separate between cases where regions do and don't neighbor reflector regions - mask = np.full(self._coremap.shape, False, dtype=bool) - mask[:,1:,:] = True - boundary = np.where(is_accel & adj_reflector & mask) - boundary_grps = boundary + (slice(None),) - net_current = net_current_back[boundary_grps] - dtilde = self._dtilde[boundary + (slice(None), 2)] - flux = cell_flux[boundary_grps] - self._dhat2[boundary_grps + (2,)] = (net_current + dtilde * flux) / flux - - boundary = np.where(is_accel & ~adj_reflector & mask) - boundary_grps = boundary + (slice(None),) - net_current = net_current_back[boundary_grps] - dtilde = self._dtilde[boundary + (slice(None), 2)] - flux = cell_flux[boundary_grps] - flux_back = neig_flux[boundary_grps] - self._dhat2[boundary_grps + (2,)] = ((net_current - dtilde * (flux_back - flux)) - / (flux_back + flux)) - - # Logical for whether neighboring cell to the front is reflector region - adj_reflector = np.roll(self._coremap, -1, axis=1) == _CMFD_NOACCEL - # Cell flux of neighbor to front - neig_flux = np.roll(self._flux, -1, axis=1) / dxdydz - - # Define dhat at front surface for all mesh cells not on front boundary - # Separate between cases where regions do and don't neighbor reflector regions - mask = np.full(self._coremap.shape, False, dtype=bool) - mask[:,:-1,:] = True - boundary = np.where(is_accel & adj_reflector & mask) - boundary_grps = boundary + (slice(None),) - net_current = net_current_front[boundary_grps] - dtilde = self._dtilde[boundary + (slice(None), 3)] - flux = cell_flux[boundary_grps] - self._dhat2[boundary_grps + (3,)] = (net_current - dtilde * flux) / flux - - boundary = np.where(is_accel & ~adj_reflector & mask) - boundary_grps = boundary + (slice(None),) - net_current = net_current_front[boundary_grps] - dtilde = self._dtilde[boundary + (slice(None), 3)] - flux = cell_flux[boundary_grps] - flux_front = neig_flux[boundary_grps] - self._dhat2[boundary_grps + (3,)] = ((net_current + dtilde * (flux_front - flux)) - / (flux_front + flux)) - - # Logical for whether neighboring cell to the bottom is reflector region - adj_reflector = np.roll(self._coremap, 1, axis=2) == _CMFD_NOACCEL - # Cell flux of neighbor to bottom - neig_flux = np.roll(self._flux, 1, axis=2) / dxdydz - - # Define dhat at bottom surface for all mesh cells not on bottom boundary - # Separate between cases where regions do and don't neighbor reflector regions - mask = np.full(self._coremap.shape, False, dtype=bool) - mask[:,:,1:] = True - boundary = np.where(is_accel & adj_reflector & mask) - boundary_grps = boundary + (slice(None),) - net_current = net_current_bottom[boundary_grps] - dtilde = self._dtilde[boundary + (slice(None), 4)] - flux = cell_flux[boundary_grps] - self._dhat2[boundary_grps + (4,)] = (net_current + dtilde * flux) / flux - - boundary = np.where(is_accel & ~adj_reflector & mask) - boundary_grps = boundary + (slice(None),) - net_current = net_current_bottom[boundary_grps] - dtilde = self._dtilde[boundary + (slice(None), 4)] - flux = cell_flux[boundary_grps] - flux_bottom = neig_flux[boundary_grps] - self._dhat2[boundary_grps + (4,)] = ((net_current - dtilde * (flux_bottom - flux)) - / (flux_bottom + flux)) - - # Logical for whether neighboring cell to the top is reflector region - adj_reflector = np.roll(self._coremap, -1, axis=2) == _CMFD_NOACCEL - # Cell flux of neighbor to top - neig_flux = np.roll(self._flux, -1, axis=2) / dxdydz - - # Define dhat at top surface for all mesh cells not on top boundary - # Separate between cases where regions do and don't neighbor reflector regions - mask = np.full(self._coremap.shape, False, dtype=bool) - mask[:,:,:-1] = True - boundary = np.where(is_accel & adj_reflector & mask) - boundary_grps = boundary + (slice(None),) - net_current = net_current_top[boundary_grps] - dtilde = self._dtilde[boundary + (slice(None), 5)] - flux = cell_flux[boundary_grps] - self._dhat2[boundary_grps + (5,)] = (net_current - dtilde * flux) / flux - - boundary = np.where(is_accel & ~adj_reflector & mask) - boundary_grps = boundary + (slice(None),) - net_current = net_current_top[boundary_grps] - dtilde = self._dtilde[boundary + (slice(None), 5)] - flux = cell_flux[boundary_grps] - flux_top = neig_flux[boundary_grps] - self._dhat2[boundary_grps + (5,)] = ((net_current + dtilde * (flux_top - flux)) - / (flux_top + flux)) - print(np.where(self._dhat2 != self._dhat)) - - def _compute_dhat(self): """Computes the nonlinear coupling coefficient by looping over all spatial regions and energy groups @@ -4046,101 +4255,4 @@ class CMFDRun(object): else: return current[2*l+1]/current[2*l] - def _create_cmfd_tally(self): - """Creates all tallies in-memory that are used to solve CMFD problem""" - # Create Mesh object based on CMFDMesh, stored internally - cmfd_mesh = openmc.capi.Mesh() - # Store id of Mesh object - self._cmfd_mesh_id = cmfd_mesh.id - # Set dimension and parameters of Mesh object - cmfd_mesh.dimension = self._cmfd_mesh.dimension - cmfd_mesh.set_parameters(lower_left=self._cmfd_mesh.lower_left, - upper_right=self._cmfd_mesh.upper_right, - width=self._cmfd_mesh.width) - # Create Mesh Filter object, stored internally - mesh_filter = openmc.capi.MeshFilter() - # Set mesh for Mesh Filter - mesh_filter.mesh = cmfd_mesh - - # Set up energy filters, if applicable - if self._energy_filters: - # Create Energy Filter object, stored internally - energy_filter = openmc.capi.EnergyFilter() - # Set bins for Energy Filter - energy_filter.bins = self._egrid - - # Create Energy Out Filter object, stored internally - energyout_filter = openmc.capi.EnergyoutFilter() - # Set bins for Energy Filter - energyout_filter.bins = self._egrid - - # Create Mesh Surface Filter object, stored internally - meshsurface_filter = openmc.capi.MeshSurfaceFilter() - # Set mesh for Mesh Surface Filter - meshsurface_filter.mesh = cmfd_mesh - - # Create Legendre Filter object, stored internally - legendre_filter = openmc.capi.LegendreFilter() - # Set order for Legendre Filter - legendre_filter.order = 1 - - # Create CMFD tallies, stored internally - n_tallies = 4 - self._cmfd_tally_ids = [] - for i in range(n_tallies): - tally = openmc.capi.Tally() - # Set nuclide bins - tally.nuclides = ['total'] - self._cmfd_tally_ids.append(tally.id) - - # Set attributes of CMFD flux, total tally - if i == 0: - # Set filters for tally - if self._energy_filters: - tally.filters = [mesh_filter, energy_filter] - else: - tally.filters = [mesh_filter] - # Set scores, type, and estimator for tally - tally.scores = ['flux', 'total'] - tally.type = 'volume' - tally.estimator = 'analog' - - # Set attributes of CMFD neutron production tally - elif i == 1: - # Set filters for tally - if self._energy_filters: - tally.filters = [mesh_filter, energy_filter, energyout_filter] - else: - tally.filters = [mesh_filter] - # Set scores, type, and estimator for tally - tally.scores = ['nu-scatter', 'nu-fission'] - tally.type = 'volume' - tally.estimator = 'analog' - - # Set attributes of CMFD surface current tally - elif i == 2: - # Set filters for tally - if self._energy_filters: - tally.filters = [meshsurface_filter, energy_filter] - else: - tally.filters = [meshsurface_filter] - # Set scores, type, and estimator for tally - tally.scores = ['current'] - tally.type = 'mesh-surface' - tally.estimator = 'analog' - - # Set attributes of CMFD P1 scatter tally - elif i == 3: - # Set filters for tally - if self._energy_filters: - tally.filters = [mesh_filter, legendre_filter, energy_filter] - else: - tally.filters = [mesh_filter, legendre_filter] - # Set scores for tally - tally.scores = ['scatter'] - tally.type = 'volume' - tally.estimator = 'analog' - - # Set all tallies to be active from beginning - tally.active = True