diff --git a/DG_Approximation.py b/DG_Approximation.py index 42ce8e44d46deeb5c50f5cf317455a8d8260675c..7412de8f930d36df636fc8bbf4d6d79e93a524c2 100644 --- a/DG_Approximation.py +++ b/DG_Approximation.py @@ -25,6 +25,9 @@ TODO: Contemplate moving plots to pertaining files TODO: Check time efficiency of details plot TODO: Contemplate moving A and B to Vectors_of_Polynomials TODO: Combine plot for coarse and fine approximation for wavelet detectors +TODO: Remove unnecessary uses of len() -> Done +TODO: Remove unnecessary uses of int() -> Done +TODO: Rename self.num_grid_cells in Troubled_Cell_Detector to self.num_coarse_grid_cells """ import numpy as np @@ -85,10 +88,6 @@ class DGScheme(object): if not hasattr(Update_Scheme, self.update_scheme): raise ValueError('Invalid update scheme: "%s"' % self.update_scheme) - # Adjust config for detectors using wavelets - if self.detector in ['Boxplot', 'Theoretical']: - self.detector_config['polynom_degree'] = self.polynom_degree - self._reset() # Replace the string names with the actual class instances @@ -272,7 +271,7 @@ class DGScheme(object): for point in range(num_points)] for cell in range(len(projection[0]))] - return np.reshape(np.array(approx), (1, len(approx) * len(approx[0]))) + return np.reshape(np.array(approx), (1, len(approx) * num_points)) def _calculate_exact_solution(self, mesh, cell_len): grid = [] @@ -305,7 +304,7 @@ class DGScheme(object): # Calculate projection on coarse mesh output_matrix = [] - for i in range(int(len(projection[0])/2)): + for i in range(self.num_grid_cells//2): new_entry = 0.5*(projection[:, 2*i] @ basis_matrix_left + projection[:, 2*i+1] @ basis_matrix_right) output_matrix.append(new_entry) coarse_projection = np.transpose(np.array(output_matrix)) @@ -371,7 +370,7 @@ class DGScheme(object): coarse_projection = self._calculate_coarse_projection(projection) averaged_projection = [[coarse_projection[degree][cell]*self.basis[degree].subs(x, value) - for cell in range(len(coarse_projection[0])) + for cell in range(self.num_grid_cells//2) for value in [-0.5, 0.5]] for degree in range(self.polynom_degree+1)] @@ -379,7 +378,7 @@ class DGScheme(object): wavelet = AlpertsWavelet(self.polynom_degree).get_vector(xi) wavelet_projection = [[multiwavelet_coeffs[degree][cell]*wavelet[degree].subs(xi, 0.5) * value - for cell in range(len(coarse_projection[0])) + for cell in range(self.num_grid_cells//2) for value in [(-1)**(self.polynom_degree+degree+1), 1]] for degree in range(self.polynom_degree+1)] @@ -402,24 +401,10 @@ class DGScheme(object): print(avgMatrix == averaged_projection) print(testMatrix == wavelet_projection) - # uAvg = np.sum(avgMatrix, axis=0) - # uH = np.sum([newMatrix[degree] * phiVector[degree].subs(x, 0) - # for degree in range(p + 1)], axis=0) - # testD = np.sum(testMatrix, axis=0) - projected_coarse = np.sum(averaged_projection, axis=0) - # print(type(projected_coarse)) - # print(projected_coarse) projected_fine = np.sum([fine_projection[degree] * self.basis[degree].subs(x, 0) for degree in range(self.polynom_degree+1)], axis=0) - # print(type(projected_fine)) - # print(projected_fine) - # print('Proj: ', projected_coarse.shape, projected_fine.shape) projected_wavelet_coeffs = np.sum(wavelet_projection, axis=0) - # print(type(projected_wavelet_coeffs)) - # print(projected_wavelet_coeffs) - # print('Wave: ', projected_wavelet_coeffs.shape) - # print(fine_mesh.shape) plt.figure(4) plt.plot(fine_mesh, projected_fine-projected_coarse, 'm-.') @@ -428,7 +413,6 @@ class DGScheme(object): plt.xlabel('X') plt.ylabel('Detail Coefficients') plt.title('Wavelet Coefficients') - plt.show() def _plot_shock_tube(self): time_history = self.update_scheme.get_time_history() diff --git a/Troubled_Cell_Detector.py b/Troubled_Cell_Detector.py index 3ec9e0de67054257ea9a84301889d46b69875c72..c81f23c0065a8bbcf7d6328b801e414f16b6ee1e 100644 --- a/Troubled_Cell_Detector.py +++ b/Troubled_Cell_Detector.py @@ -39,6 +39,7 @@ class WaveletDetector(TroubledCellDetector): def __init__(self, config): # Unpack necessary configurations self.polynom_degree = config.pop('polynom_degree') + self.num_grid_cells = config.pop('num_grid_cells')//2 # Set fixed basis and wavelet vectors self.basis = OrthonormalLegendre(self.polynom_degree).get_vector(x) @@ -63,17 +64,11 @@ class WaveletDetector(TroubledCellDetector): def _calculate_wavelet_coeffs(self, projection): output_matrix = [] - for i in range(int(len(projection[0])/2)): + for i in range(self.num_grid_cells): new_entry = 0.5*(projection[:, 2*i] @ self.M1 + projection[:, 2*i+1] @ self.M2) output_matrix.append(new_entry) return np.transpose(np.array(output_matrix)) - # def get_wavelet_coeffs(self, projection=None): - # if projection is None: - # return self.multiwavelet_coeffs - # else: - # return self._calculate_wavelet_coeffs(projection) - def get_cells(self, projection): multiwavelet_coeffs = self._calculate_wavelet_coeffs(projection[:, 1: -1]) return self._get_cells(multiwavelet_coeffs, projection) @@ -90,30 +85,26 @@ class Boxplot(WaveletDetector): self.function_name = 'Boxplot' # Unpack necessary configurations - self.num_grid_cells = config.pop('num_grid_cells') self.fold_len = config.pop('fold_len', 16) self.whisker_len = config.pop('whisker_len', 3) # Pop unnecessary parameter - config.pop('polynom_degree') config.pop('cell_len') def _get_cells(self, multiwavelet_coeffs, projection): - self.num_grid_cells = len(multiwavelet_coeffs[0]) indexed_coeffs = [[multiwavelet_coeffs[0, i], i]for i in range(self.num_grid_cells)] if self.num_grid_cells < self.fold_len: self.fold_len = self.num_grid_cells - num_folds = int(self.num_grid_cells/self.fold_len) + num_folds = self.num_grid_cells//self.fold_len troubled_cells = [] for fold in range(num_folds): sorted_fold = sorted(indexed_coeffs[fold * self.fold_len:(fold+1) * self.fold_len]) - # fold_len/4.0 had comment: int((P+3)/2)/2.0 - balance_factor = self.fold_len/4.0 % 1 # former: int(...) - boundary_index = int(self.fold_len/4.0 - balance_factor) + boundary_index = self.fold_len//4 + balance_factor = self.fold_len/4.0 - boundary_index first_quartile = (1-balance_factor) * sorted_fold[boundary_index-1][0] \ + balance_factor * sorted_fold[boundary_index][0] @@ -148,7 +139,6 @@ class Theoretical(WaveletDetector): self.function_name = 'Theoretical' # Unpack necessary configurations - self.num_grid_cells = config.pop('num_grid_cells') self.cell_len = config.pop('cell_len') self.cutoff_factor = config.pop('cutoff_factor', np.sqrt(2) * self.cell_len) # comment to line above: or 2 or 3 @@ -161,7 +151,7 @@ class Theoretical(WaveletDetector): troubled_cells = [] self.max_avg = max(1, max(abs(projection[0][degree]) for degree in range(self.polynom_degree+1))) - for cell in range(int(self.num_grid_cells/2)): + for cell in range(self.num_grid_cells): if self._is_troubled_cell(cell): troubled_cells.append(cell) diff --git a/Vectors_of_Polynomials.py b/Vectors_of_Polynomials.py index 7ce26f8ba4707ea63d025d5af0e0d70d1828cb95..07b597fbe3d53066a360fa9c4ef7037407313185 100644 --- a/Vectors_of_Polynomials.py +++ b/Vectors_of_Polynomials.py @@ -27,8 +27,7 @@ class Legendre(Vector): if degree == 1: vector.append(eval_point) else: - poly = (2.0*degree - 1)/degree * eval_point * vector[len(vector)-1] \ - - (degree-1)/degree * vector[len(vector)-2] + poly = (2.0*degree - 1)/degree * eval_point * vector[-1] - (degree-1)/degree * vector[-2] vector.append(poly) return vector @@ -36,7 +35,7 @@ class Legendre(Vector): class OrthonormalLegendre(Legendre): def get_vector(self, eval_point): leg_vector = self._calculate_legendre_vector(eval_point) - return [leg_vector[degree] * np.sqrt(degree+0.5) for degree in range(len(leg_vector))] + return [leg_vector[degree] * np.sqrt(degree+0.5) for degree in range(self.polynom_degree+1)] class AlpertsWavelet(Vector):