diff --git a/Limiter.py b/Limiter.py index 892a4fe4f9e4ca5741c9ebaa3032ab053368a8c9..22e0189f31484e23be8064c6db876fe01be2babd 100644 --- a/Limiter.py +++ b/Limiter.py @@ -45,6 +45,13 @@ class MinMod(Limiter): adapted_projection[i] = 0 return adapted_projection + def _determine_modification(self, projection, cell, cell_slope): + forward_slope = (projection[0][cell+1] - projection[0][cell]) * (0.5**0.5) + backward_slope = (projection[0][cell] - projection[0][cell-1]) * (0.5**0.5) + + return (forward_slope <= 0) & (backward_slope <= 0) & (cell_slope <= 0) \ + | (forward_slope >= 0) & (backward_slope >= 0) & (cell_slope >= 0) + @staticmethod def _set_cell_slope(projection, cell): slope = [] @@ -54,13 +61,6 @@ class MinMod(Limiter): slope.append(new_entry) return slope[cell] - def _determine_modification(self, projection, cell, cell_slope): - forward_slope = (projection[0][cell+1] - projection[0][cell]) * (0.5**0.5) - backward_slope = (projection[0][cell] - projection[0][cell-1]) * (0.5**0.5) - - return (forward_slope <= 0) & (backward_slope <= 0) & (cell_slope <= 0) \ - | (forward_slope >= 0) & (backward_slope >= 0) & (cell_slope >= 0) - class ModifiedMinMod(MinMod): def _reset(self, config): diff --git a/Quadrature.py b/Quadrature.py index fc3babfac75b4bfdd2f64c084113206b6801a83b..600abfd4cc48bf64fc4e047eddf1db177fdf230a 100644 --- a/Quadrature.py +++ b/Quadrature.py @@ -11,22 +11,22 @@ class Quadrature(object): self._reset(config) def _reset(self, config): + self._num_eval_points = None self._eval_points = None self._weights = None - self._num_eval_points = None def get_name(self): return self.__class__.__name__ + def get_num_points(self): + return self._num_eval_points + def get_eval_points(self): return self._eval_points def get_weights(self): return self._weights - def get_num_points(self): - return self._num_eval_points - class Gauss(Quadrature): def _reset(self, config): diff --git a/Troubled_Cell_Detector.py b/Troubled_Cell_Detector.py index 8adf214b07e8724ea06e62d66d959311da19183c..f21455728b63a4a489d1fc100420c4e3bf74f79c 100644 --- a/Troubled_Cell_Detector.py +++ b/Troubled_Cell_Detector.py @@ -51,6 +51,17 @@ class TroubledCellDetector(object): print("N =", self._num_grid_cells) print("maximum error =", max_error) + def _plot_shock_tube(self, troubled_cell_history, time_history): + plt.figure(6) + for pos in range(len(time_history)): + current_cells = troubled_cell_history[pos] + for cell in current_cells: + plt.plot(cell, time_history[pos], 'k.') + plt.xlim((0, self._num_grid_cells//2)) + plt.xlabel('Cell') + plt.ylabel('Time') + plt.title('Shock Tubes') + def _plot_mesh(self, projection, color_exact, color_approx): grid, exact = self._calculate_exact_solution(self._mesh[2:-2], self._cell_len) approx = self._calculate_approximate_solution(projection[:, 1:-1]) @@ -88,17 +99,6 @@ class TroubledCellDetector(object): plt.ylabel('u(x,t)-uh(x,t)') plt.title('Errors') - def _plot_shock_tube(self, troubled_cell_history, time_history): - plt.figure(6) - for pos in range(len(time_history)): - current_cells = troubled_cell_history[pos] - for cell in current_cells: - plt.plot(cell, time_history[pos], 'k.') - plt.xlim((0, self._num_grid_cells//2)) - plt.xlabel('Cell') - plt.ylabel('Time') - plt.title('Shock Tubes') - def _calculate_exact_solution(self, mesh, cell_len): grid = [] exact = [] @@ -170,9 +170,6 @@ class WaveletDetector(TroubledCellDetector): multiwavelet_coeffs = self._calculate_wavelet_coeffs(projection[:, 1: -1]) return self._get_cells(multiwavelet_coeffs, projection) - def _get_cells(self, multiwavelet_coeffs, projection): - return [] - def _calculate_wavelet_coeffs(self, projection): output_matrix = [] for i in range(self._num_coarse_grid_cells): @@ -180,37 +177,13 @@ class WaveletDetector(TroubledCellDetector): output_matrix.append(new_entry) return np.transpose(np.array(output_matrix)) + def _get_cells(self, multiwavelet_coeffs, projection): + return [] + def plot_results(self, projection, troubled_cell_history, time_history, color_exact, color_approx): self._plot_details(projection) super().plot_results(projection, troubled_cell_history, time_history, color_exact, color_approx) - def _plot_mesh(self, projection, color_exact, color_approx): - grid, exact = self._calculate_exact_solution(self._mesh[2:-2], self._cell_len) - approx = self._calculate_approximate_solution(projection[:, 1:-1]) - - pointwise_error = np.abs(exact-approx) - max_error = np.max(pointwise_error) - - self._plot_coarse_mesh(projection, color_exact, color_approx) - self._plot_solution_and_approx(grid, exact, approx, 'k-.', 'b-.') - plt.legend(['Exact (Coarse)', 'Approx (Coarse)', 'Exact (Fine)', 'Approx (Fine)']) - self._plot_semilog_error(grid, pointwise_error) - self._plot_error(grid, exact, approx) - - return max_error - - def _plot_coarse_mesh(self, projection, color_exact, color_approx): - coarse_cell_len = 2*self._cell_len - coarse_mesh = np.arange(self._left_bound - (0.5*coarse_cell_len), self._right_bound + (1.5*coarse_cell_len), - coarse_cell_len) - - coarse_projection = self._calculate_coarse_projection(projection) - - # Plot exact and approximate solutions for coarse mesh - grid, exact = self._calculate_exact_solution(coarse_mesh[1:-1], coarse_cell_len) - approx = self._calculate_approximate_solution(coarse_projection) - self._plot_solution_and_approx(grid, exact, approx, color_exact, color_approx) - def _plot_details(self, projection): fine_mesh = self._mesh[2:-2] @@ -274,8 +247,8 @@ class WaveletDetector(TroubledCellDetector): plt.title('Wavelet Coefficients') def _calculate_coarse_projection(self, projection): - basis_matrix_left = self._set_basis_matrix(xi, 0.5*(xi-1)) - basis_matrix_right = self._set_basis_matrix(xi, 0.5*(xi+1)) + basis_matrix_left = self._set_basis_matrix(xi, 0.5 * (xi - 1)) + basis_matrix_right = self._set_basis_matrix(xi, 0.5 * (xi + 1)) # Remove ghost cells projection = projection[:, 1:-1] @@ -283,7 +256,8 @@ class WaveletDetector(TroubledCellDetector): # Calculate projection on coarse mesh output_matrix = [] for i in range(self._num_coarse_grid_cells): - new_entry = 0.5*(projection[:, 2*i] @ basis_matrix_left + projection[:, 2*i+1] @ basis_matrix_right) + 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)) @@ -291,15 +265,42 @@ class WaveletDetector(TroubledCellDetector): def _set_basis_matrix(self, first_param, second_param): matrix = [] - for i in range(self._polynom_degree+1): + for i in range(self._polynom_degree + 1): row = [] - for j in range(self._polynom_degree+1): + for j in range(self._polynom_degree + 1): entry = integrate(self._basis[i].subs(x, first_param) * self._basis[j].subs(x, second_param), (xi, -1, 1)) row.append(np.float64(entry)) matrix.append(row) return matrix + def _plot_mesh(self, projection, color_exact, color_approx): + grid, exact = self._calculate_exact_solution(self._mesh[2:-2], self._cell_len) + approx = self._calculate_approximate_solution(projection[:, 1:-1]) + + pointwise_error = np.abs(exact-approx) + max_error = np.max(pointwise_error) + + self._plot_coarse_mesh(projection, color_exact, color_approx) + self._plot_solution_and_approx(grid, exact, approx, 'k-.', 'b-.') + plt.legend(['Exact (Coarse)', 'Approx (Coarse)', 'Exact (Fine)', 'Approx (Fine)']) + self._plot_semilog_error(grid, pointwise_error) + self._plot_error(grid, exact, approx) + + return max_error + + def _plot_coarse_mesh(self, projection, color_exact, color_approx): + coarse_cell_len = 2*self._cell_len + coarse_mesh = np.arange(self._left_bound - (0.5*coarse_cell_len), self._right_bound + (1.5*coarse_cell_len), + coarse_cell_len) + + coarse_projection = self._calculate_coarse_projection(projection) + + # Plot exact and approximate solutions for coarse mesh + grid, exact = self._calculate_exact_solution(coarse_mesh[1:-1], coarse_cell_len) + approx = self._calculate_approximate_solution(coarse_projection) + self._plot_solution_and_approx(grid, exact, approx, color_exact, color_approx) + class Boxplot(WaveletDetector): def _reset(self, config): diff --git a/Update_Scheme.py b/Update_Scheme.py index 2f10cd2c832bfcb1550706c6a0f00c922c6437d1..42464f4266def1c6c6037090b0e3e61fc3273ac9 100644 --- a/Update_Scheme.py +++ b/Update_Scheme.py @@ -26,17 +26,6 @@ class UpdateScheme(object): self._reset() - def get_name(self): - return self.__class__.__name__ - - def step(self, projection, cfl_number): - current_projection, troubled_cells = self._apply_stability_method(projection, cfl_number) - - return current_projection, troubled_cells - - def _apply_stability_method(self, projection, cfl_number): - return projection, [] - def _reset(self): # Set matrix A matrix = [] @@ -60,6 +49,17 @@ class UpdateScheme(object): matrix.append(new_row) self._B = np.array(matrix) # former: inv_mass @ np.array(matrix) + def get_name(self): + return self.__class__.__name__ + + def step(self, projection, cfl_number): + current_projection, troubled_cells = self._apply_stability_method(projection, cfl_number) + + return current_projection, troubled_cells + + def _apply_stability_method(self, projection, cfl_number): + return projection, [] + def _apply_limiter(self, current_projection): troubled_cells = self._detector.get_cells(current_projection) @@ -94,6 +94,18 @@ class SSPRK3(UpdateScheme): return current_projection, troubled_cells + def _apply_first_step(self, original_projection, cfl_number): + right_hand_side = self._update_right_hand_side(original_projection) + return original_projection + (cfl_number*right_hand_side) + + def _apply_second_step(self, original_projection, current_projection, cfl_number): + right_hand_side = self._update_right_hand_side(current_projection) + return 1/4 * (3*original_projection + (current_projection + cfl_number*right_hand_side)) + + def _apply_third_step(self, original_projection, current_projection, cfl_number): + right_hand_side = self._update_right_hand_side(current_projection) + return 1/3 * (original_projection + 2*(current_projection + cfl_number*right_hand_side)) + def _update_right_hand_side(self, current_projection): # Initialize vector and set first entry to accommodate for ghost cell right_hand_side = [0] @@ -107,15 +119,3 @@ class SSPRK3(UpdateScheme): right_hand_side.append(right_hand_side[1]) return np.transpose(right_hand_side) - - def _apply_first_step(self, original_projection, cfl_number): - right_hand_side = self._update_right_hand_side(original_projection) - return original_projection + (cfl_number*right_hand_side) - - def _apply_second_step(self, original_projection, current_projection, cfl_number): - right_hand_side = self._update_right_hand_side(current_projection) - return 1/4 * (3*original_projection + (current_projection + cfl_number*right_hand_side)) - - def _apply_third_step(self, original_projection, current_projection, cfl_number): - right_hand_side = self._update_right_hand_side(current_projection) - return 1/3 * (original_projection + 2*(current_projection + cfl_number*right_hand_side))