From 2313d35d165c9f8ceaa5d04ed7a2e6770fde4eb3 Mon Sep 17 00:00:00 2001 From: lakue103 <laura.kuehle@uni-duesseldorf.de> Date: Wed, 2 Sep 2020 13:13:37 +0200 Subject: [PATCH] Moved everything pertaining wavelets to Troubled_Cell_Detector --- DG_Approximation.py | 34 +++++++++++------- Main.py | 3 ++ Troubled_Cell_Detector.py | 75 +++++++++++++++++++++++++++++++++++---- Update_Scheme.py | 59 +++--------------------------- 4 files changed, 98 insertions(+), 73 deletions(-) diff --git a/DG_Approximation.py b/DG_Approximation.py index e94626c..1a8f52f 100644 --- a/DG_Approximation.py +++ b/DG_Approximation.py @@ -25,6 +25,9 @@ TODO: Write documentation for all methods TODO: Add a verbose option TODO: Check whether consistency is given/possible for each class instance +TODO: Make projection local variable +TODO: Check whether current and original projection in Update_Scheme are identical + """ import numpy as np from sympy import Symbol, integrate @@ -90,6 +93,10 @@ class DGScheme(object): 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 @@ -135,16 +142,19 @@ class DGScheme(object): self._plot_shock_tube() self._plot_coarse_mesh('k-', 'y') - approx = self._plot_fine_mesh('k-.', 'b') - troubled_cells = self.update_scheme.get_troubled_cell_history() - wavelet = self.update_scheme.get_wavelet_coeffs( - self.projection[:, 1: -1]) + approx = self._plot_fine_mesh('k-.', 'b') - return approx + # What is that??? What is it for? - return self.projection[:, 1:-1], self.mesh[2:-2], troubled_cells, \ - wavelet + # troubled_cells = self.update_scheme.get_troubled_cell_history() + # wavelet = self.update_scheme.get_wavelet_coeffs( + # self.projection[:, 1: -1]) + # + # return approx + # + # return self.projection[:, 1:-1], self.mesh[2:-2], troubled_cells, \ + # wavelet def save_plots(self): name = self.init_cond.get_name() + '__' \ @@ -245,7 +255,7 @@ class DGScheme(object): return matrix def _do_initial_projection(self): - # Initialze matrix and set first entry to accommodate for ghost cell + # Initialize matrix and set first entry to accommodate for ghost cell output_matrix = [0] for cell in range(self.num_grid_cells): @@ -274,7 +284,7 @@ class DGScheme(object): print(np.array(output_matrix).shape) self.projection = np.transpose(np.array(output_matrix)) - def _calculate_approximative_solution(self, projection): + def _calculate_approximate_solution(self, projection): points = self.quadrature.get_eval_points() num_points = self.quadrature.get_num_points() basis = self.basis @@ -357,10 +367,10 @@ class DGScheme(object): print(coarse_projection == coarse_projection1) - # Plot exact and approximative solutions for coarse mesh + # Plot exact and approximate solutions for coarse mesh grid, exact = self._calculate_exact_solution(coarse_mesh[1:-1], coarse_cell_len) - approx = self._calculate_approximative_solution(coarse_projection) + approx = self._calculate_approximate_solution(coarse_projection) self._plot_solution_and_approx(grid, exact, approx, color_exact, color_approx) # plt.legend(['Exact-Coarse', 'Approx-Coarse']) @@ -368,7 +378,7 @@ class DGScheme(object): def _plot_fine_mesh(self, color_exact, color_approx): grid, exact = self._calculate_exact_solution(self.mesh[2:-2], self.cell_len) - approx = self._calculate_approximative_solution( + approx = self._calculate_approximate_solution( self.projection[:, 1:-1]) pointwise_error = np.abs(exact-approx) diff --git a/Main.py b/Main.py index 098ba84..8796233 100644 --- a/Main.py +++ b/Main.py @@ -64,3 +64,6 @@ dg_scheme.save_plots() # print(dg_scheme.limiter) # # print(dg_scheme.projection) # ============================================================================= + +# if __name__ == '__main__': +# pass diff --git a/Troubled_Cell_Detector.py b/Troubled_Cell_Detector.py index 4013f95..7b0a4d8 100644 --- a/Troubled_Cell_Detector.py +++ b/Troubled_Cell_Detector.py @@ -7,16 +7,23 @@ TODO: Check cutoff_factor/whether _is_troubled_cell works correctly \ """ import numpy as np +from sympy import Symbol, integrate + +from Vectors_of_Polynomials import OrthonormalLegendre, AlpertsWavelet + +x = Symbol('x') +xi = Symbol('z') class TroubledCellDetector(object): def __init__(self, config): + self.function_name = 'None' pass def get_name(self): return self.function_name - def get_cells(self, multiwavelet_coeffs, projection): + def get_cells(self, projection): pass @@ -26,12 +33,65 @@ class NoDetection(TroubledCellDetector): self.function_name = 'NoDetection' pass - def get_cells(self, multiwavelet_coeffs, projection): + def get_cells(self, projection): + return [] + + +class WaveletDetector(TroubledCellDetector): + def __init__(self, config): + # Unpack necessary configurations + self.polynom_degree = config.pop('polynom_degree') + + # Set fixed basis and wavelet vectors + self.basis = OrthonormalLegendre( + self.polynom_degree).get_vector(x) + self.wavelet = AlpertsWavelet(self.polynom_degree).get_vector(x) + + # Set additional necessary parameter + self.M1 = self._set_multiwavelet_matrix(xi, -0.5*(xi-1), True) + self.M2 = self._set_multiwavelet_matrix(xi, 0.5*(xi+1), False) + + def _set_multiwavelet_matrix(self, first_param, second_param, is_M1): + matrix = [] + for i in range(self.polynom_degree+1): + row = [] + for j in range(self.polynom_degree+1): + entry = integrate(self.basis[i].subs(x, first_param) + * self.wavelet[j].subs(x, second_param), + (xi, -1, 1)) + if is_M1: + entry = entry * (-1)**(j + self.polynom_degree + 1) + row.append(np.float64(entry)) + matrix.append(row) + return matrix + + def _calculate_wavelet_coeffs(self, projection): + transposed_vector = np.transpose(projection) + output_matrix = [] + for i in range(int(len(projection[0])/2)): + new_entry = 0.5*(transposed_vector[2*i] @ self.M1 + + transposed_vector[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) + + def _get_cells(self, multiwavelet_coeffs, projection): return [] -class Boxplot(TroubledCellDetector): +class Boxplot(WaveletDetector): def __init__(self, config): + super().__init__(config) + # Set name of function self.function_name = 'Boxplot' @@ -44,7 +104,7 @@ class Boxplot(TroubledCellDetector): config.pop('polynom_degree') config.pop('cell_len') - def get_cells(self, multiwavelet_coeffs, projection): + 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)] @@ -92,13 +152,14 @@ class Boxplot(TroubledCellDetector): return sorted(troubled_cells) -class Theoretical(TroubledCellDetector): +class Theoretical(WaveletDetector): def __init__(self, config): + super().__init__(config) + # Set name of function self.function_name = 'Theoretical' # Unpack necessary configurations - self.polynom_degree = config.pop('polynom_degree') self.num_grid_cells = config.pop('num_grid_cells') self.cell_len = config.pop('cell_len') self.cutoff_factor = config.pop('cutoff_factor', @@ -109,7 +170,7 @@ class Theoretical(TroubledCellDetector): self.projection = [] self.max_avg = None - def get_cells(self, multiwavelet_coeffs, projection): + def _get_cells(self, multiwavelet_coeffs, projection): self.multiwavelet_coeffs = multiwavelet_coeffs self.projection = projection troubled_cells = [] diff --git a/Update_Scheme.py b/Update_Scheme.py index 097389a..1351307 100644 --- a/Update_Scheme.py +++ b/Update_Scheme.py @@ -13,12 +13,6 @@ TODO: Find better names for A, B, M1, and M2 """ import numpy as np -from sympy import Symbol, integrate - -from Vectors_of_Polynomials import OrthonormalLegendre, AlpertsWavelet - -x = Symbol('x') -xi = Symbol('z') class UpdateScheme(object): @@ -50,12 +44,6 @@ class UpdateScheme(object): def get_time_history(self): return self.time_history - def get_wavelet_coeffs(self, projection=None): - if projection is None: - return self.multiwavelet_coeffs - else: - return self._calculate_wavelet_coeffs(projection) - def step(self, projection, cfl_number, time): self.original_projection = projection self.current_projection = projection @@ -76,17 +64,10 @@ class UpdateScheme(object): pass def _reset(self): - # Set fixed basis and wavelet vectors - self.basis = OrthonormalLegendre( - self.polynom_degree).get_vector(x) - self.wavelet = AlpertsWavelet(self.polynom_degree).get_vector(x) - # Set additional necessary fixed instance variables self.name = 'None' self.interval_len = self.right_bound-self.left_bound self.cell_len = self.interval_len / self.num_grid_cells - self.M1 = self._set_multiwavelet_matrix(xi, -0.5*(xi-1), True) - self.M2 = self._set_multiwavelet_matrix(xi, 0.5*(xi+1), False) # Set matrix A matrix = [] @@ -114,7 +95,6 @@ class UpdateScheme(object): self.original_projection = [] self.current_projection = [] self.right_hand_side = [] - self.multiwavelet_coeffs = [] self.troubled_cells = [] self.troubled_cell_history = [] self.time_history = [] @@ -122,37 +102,8 @@ class UpdateScheme(object): self.time = 0 self.iteration = 0 - def _set_multiwavelet_matrix(self, first_param, second_param, is_M1): - matrix = [] - for i in range(self.polynom_degree+1): - row = [] - for j in range(self.polynom_degree+1): - entry = integrate(self.basis[i].subs(x, first_param) - * self.wavelet[j].subs(x, second_param), - (xi, -1, 1)) - if is_M1: - entry = entry * (-1)**(j + self.polynom_degree + 1) - row.append(np.float64(entry)) - matrix.append(row) - return matrix - - def _update_wavelet_coeffs(self): - projection = self.current_projection[:, 1: -1] - self.multiwavelet_coeffs = self._calculate_wavelet_coeffs(projection) - - def _calculate_wavelet_coeffs(self, projection): - transposed_vector = np.transpose(projection) - output_matrix = [] - for i in range(int(len(projection[0])/2)): - new_entry = 0.5*(transposed_vector[2*i] @ self.M1 - + transposed_vector[2*i+1] @ self.M2) - output_matrix.append(new_entry) - return np.transpose(np.array(output_matrix)) - def _apply_limiter(self): - self._update_wavelet_coeffs() - self.troubled_cells = self.detector.get_cells(self.multiwavelet_coeffs, - self.current_projection) + self.troubled_cells = self.detector.get_cells(self.current_projection) new_projection = self.current_projection.copy() for cell in self.troubled_cells: @@ -174,13 +125,13 @@ class SSPRK3(UpdateScheme): def __init__(self, detector, limiter, init_cond, mesh, wave_speed, polynom_degree, num_grid_cells, final_time, history_threshold, left_bound, right_bound): - # Set name of update scheme - self.name = 'SSPRK3' - super().__init__(detector, limiter, init_cond, mesh, wave_speed, polynom_degree, num_grid_cells, final_time, history_threshold, left_bound, right_bound) + # Set name of update scheme + self.name = 'SSPRK3' + # Override method of superclass def _apply_stability_method(self): self._apply_first_step() @@ -199,7 +150,7 @@ class SSPRK3(UpdateScheme): # Transpose projection for easier calculation transposed_projection = np.transpose(self.current_projection) - # Initialze vector and set first entry to accommodate for ghost cell + # Initialize vector and set first entry to accommodate for ghost cell right_hand_side = [0] for j in range(self.num_grid_cells): -- GitLab