From d57592a0314ce910f832a02e0e9d47cacf47f87f Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Wed, 31 Oct 2018 01:18:24 -0400 Subject: [PATCH] Initial attempt at making vectorized functions more readable --- openmc/cmfd.py | 865 ++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 853 insertions(+), 12 deletions(-) diff --git a/openmc/cmfd.py b/openmc/cmfd.py index 38f67ab5a..9a82d0fd8 100644 --- a/openmc/cmfd.py +++ b/openmc/cmfd.py @@ -1302,12 +1302,14 @@ class CMFDRun(object): # Calculate dtilde if vectorized: self._compute_dtilde_vectorized() + self._compute_dtilde_vectorized_updated() else: self._compute_dtilde() # Calculate dhat if vectorized: self._compute_dhat_vectorized() + self._compute_dhat_vectorized_updated() else: self._compute_dhat() @@ -1464,12 +1466,15 @@ class CMFDRun(object): # Initialize CMFD flux map that maps phi to actual spatial and # group indices of problem cmfd_flux = np.zeros((nx, ny, nz, ng)) + cmfd_flux2 = np.zeros((nx, ny, nz, ng)) # Loop over all groups and set CMFD flux based on indices of # coremap and values of phi for g in range(ng): 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(np.where(cmfd_flux != cmfd_flux2)) # Compute fission source cmfd_src = np.sum(self._nfissxs[:,:,:,:,:] * \ @@ -1594,6 +1599,16 @@ class CMFDRun(object): energy_bins = np.where(source_energies < energy[0], ng - 1, np.where(source_energies > energy[-1], 0, \ ng - np.digitize(source_energies, energy))) + # Determine which energy bin each particle's energy belongs to + # Separate into cases bases on where source energies lies on egrid + energy_bins2 = np.zeros(len(source_energies), dtype=int) + idx = np.where(source_energies < energy[0]) + energy_bins2[idx] = ng - 1 + idx = np.where(source_energies > energy[-1]) + 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(np.where(energy_bins != energy_bins2)) # Determine weight factor of each particle based on its mesh index # and energy bin and updates its weight @@ -1644,6 +1659,17 @@ class CMFDRun(object): np.where(source_energies > energy[-1], ng - 1, \ np.digitize(source_energies, energy) - 1)) + # Determine which energy bin each particle's energy belongs to + # Separate into cases bases on where source energies lies on egrid + energy_bins2 = np.zeros(len(source_energies), dtype=int) + idx = np.where(source_energies < energy[0]) + energy_bins2[idx] = 0 + idx = np.where(source_energies > energy[-1]) + 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(np.where(energy_bins != energy_bins2)) + # Determine all unique combinations of mesh bin and energy bin, and # count number of particles that belong to these combinations idx, counts = np.unique(np.array([mesh_bins, energy_bins]), axis=1, @@ -1686,7 +1712,9 @@ class CMFDRun(object): # Build loss and production matrices if vectorized: loss = self._build_loss_matrix_vectorized(adjoint) + loss2 = self._build_loss_matrix_vectorized_updated(adjoint,loss) prod = self._build_prod_matrix_vectorized(adjoint) + prod2 = self._build_prod_matrix_vectorized_updated(adjoint,prod) else: loss = self._build_loss_matrix(adjoint) prod = self._build_prod_matrix(adjoint) @@ -2048,6 +2076,215 @@ 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 @@ -2155,6 +2392,39 @@ 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 @@ -2866,6 +3136,328 @@ 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 @@ -2982,27 +3574,27 @@ class CMFDRun(object): # Extract indices of coremap that are accelerated is_accel = self._coremap != _CMFD_NOACCEL - # Define dtilde at left surface for all mesh cells on left boundary + # Define dhat at left surface for all mesh cells on left boundary self._dhat[0,:,:,:,0] = np.where(is_accel[0,:,:,np.newaxis], (net_current_left[0,:,:,:] + self._dtilde[0,:,:,:,0] * \ cell_flux[0,:,:,:]) / cell_flux[0,:,:,:], 0) - # Define dtilde at right surface for all mesh cells on right boundary + # Define dhat at right surface for all mesh cells on right boundary self._dhat[-1,:,:,:,1] = np.where(is_accel[-1,:,:,np.newaxis], (net_current_right[-1,:,:,:] - self._dtilde[-1,:,:,:,1] * \ cell_flux[-1,:,:,:]) / cell_flux[-1,:,:,:], 0) - # Define dtilde at back surface for all mesh cells on back boundary + # Define dhat at back surface for all mesh cells on back boundary self._dhat[:,0,:,:,2] = np.where(is_accel[:,0,:,np.newaxis], (net_current_back[:,0,:,:] + self._dtilde[:,0,:,:,2] * \ cell_flux[:,0,:,:]) / cell_flux[:,0,:,:], 0) - # Define dtilde at front surface for all mesh cells on front boundary + # Define dhat at front surface for all mesh cells on front boundary self._dhat[:,-1,:,:,3] = np.where(is_accel[:,-1,:,np.newaxis], (net_current_front[:,-1,:,:] - self._dtilde[:,-1,:,:,3] * \ cell_flux[:,-1,:,:]) / cell_flux[:,-1,:,:], 0) - # Define dtilde at bottom surface for all mesh cells on bottom boundary + # Define dhat at bottom surface for all mesh cells on bottom boundary self._dhat[:,:,0,:,4] = np.where(is_accel[:,:,0,np.newaxis], (net_current_bottom[:,:,0,:] + self._dtilde[:,:,0,:,4] * \ cell_flux[:,:,0,:]) / cell_flux[:,:,0,:], 0) - # Define dtilde at top surface for all mesh cells on top boundary + # Define dhat at top surface for all mesh cells on top boundary self._dhat[:,:,-1,:,5] = np.where(is_accel[:,:,-1,np.newaxis], (net_current_top[:,:,-1,:] - self._dtilde[:,:,-1,:,5] * \ cell_flux[:,:,-1,:]) / cell_flux[:,:,-1,:], 0) @@ -3012,7 +3604,7 @@ class CMFDRun(object): # Cell flux of neighbor to left neig_flux = np.roll(self._flux, 1, axis=0) / \ np.prod(self._hxyz, axis=3)[:,:,:,np.newaxis] - # Define dtilde at left surface for all mesh cells not on left boundary + # Define dhat at left surface for all mesh cells not on left boundary self._dhat[1:,:,:,:,0] = np.where(is_accel[1:,:,:,np.newaxis], \ np.where(adj_reflector[1:,:,:,np.newaxis], (net_current_left[1:,:,:,:] + self._dtilde[1:,:,:,:,0] * \ @@ -3026,7 +3618,7 @@ class CMFDRun(object): # Cell flux of neighbor to right neig_flux = np.roll(self._flux, -1, axis=0) / \ np.prod(self._hxyz, axis=3)[:,:,:,np.newaxis] - # Define dtilde at right surface for all mesh cells not on right boundary + # Define dhat at right surface for all mesh cells not on right boundary self._dhat[:-1,:,:,:,1] = np.where(is_accel[:-1,:,:,np.newaxis], \ np.where(adj_reflector[:-1,:,:,np.newaxis], (net_current_right[:-1,:,:,:] - self._dtilde[:-1,:,:,:,1] * \ @@ -3040,7 +3632,7 @@ class CMFDRun(object): # Cell flux of neighbor to back neig_flux = np.roll(self._flux, 1, axis=1) / \ np.prod(self._hxyz, axis=3)[:,:,:,np.newaxis] - # Define dtilde at back surface for all mesh cells not on back boundary + # Define dhat at back surface for all mesh cells not on back boundary self._dhat[:,1:,:,:,2] = np.where(is_accel[:,1:,:,np.newaxis], \ np.where(adj_reflector[:,1:,:,np.newaxis], (net_current_back[:,1:,:,:] + self._dtilde[:,1:,:,:,2] * \ @@ -3054,7 +3646,7 @@ class CMFDRun(object): # Cell flux of neighbor to front neig_flux = np.roll(self._flux, -1, axis=1) / \ np.prod(self._hxyz, axis=3)[:,:,:,np.newaxis] - # Define dtilde at front surface for all mesh cells not on front boundary + # Define dhat at front surface for all mesh cells not on front boundary self._dhat[:,:-1,:,:,3] = np.where(is_accel[:,:-1,:,np.newaxis], \ np.where(adj_reflector[:,:-1,:,np.newaxis], (net_current_front[:,:-1,:,:] - self._dtilde[:,:-1,:,:,3] * \ @@ -3068,7 +3660,7 @@ class CMFDRun(object): # Cell flux of neighbor to bottom neig_flux = np.roll(self._flux, 1, axis=2) / \ np.prod(self._hxyz, axis=3)[:,:,:,np.newaxis] - # Define dtilde at bottom surface for all mesh cells not on bottom boundary + # Define dhat at bottom surface for all mesh cells not on bottom boundary self._dhat[:,:,1:,:,4] = np.where(is_accel[:,:,1:,np.newaxis], \ np.where(adj_reflector[:,:,1:,np.newaxis], (net_current_bottom[:,:,1:,:] + self._dtilde[:,:,1:,:,4] * \ @@ -3082,7 +3674,7 @@ class CMFDRun(object): # Cell flux of neighbor to top neig_flux = np.roll(self._flux, -1, axis=2) / \ np.prod(self._hxyz, axis=3)[:,:,:,np.newaxis] - # Define dtilde at top surface for all mesh cells not on top boundary + # Define dhat at top surface for all mesh cells not on top boundary self._dhat[:,:,:-1,:,5] = np.where(is_accel[:,:,:-1,np.newaxis], \ np.where(adj_reflector[:,:,:-1,np.newaxis], (net_current_top[:,:,:-1,:] - self._dtilde[:,:,:-1,:,5] * \ @@ -3091,6 +3683,255 @@ 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