diff --git a/Basis_Function.py b/Basis_Function.py index 63da57912fabd0f19c55fd66c9cd743ae8648697..4184dd0a708470a4af5c6c44dcc51782219ce55b 100644 --- a/Basis_Function.py +++ b/Basis_Function.py @@ -21,6 +21,8 @@ class Basis: Attributes ---------- + polynomial degree : int + Polynomial degree. basis : ndarray Array of basis. wavelet : ndarray @@ -47,6 +49,11 @@ class Basis: """ self._polynomial_degree = polynomial_degree + @property + def polynomial_degree(self): + """Return polynomial degree.""" + return self._polynomial_degree + @property @cache def basis(self): diff --git a/DG_Approximation.py b/DG_Approximation.py index ef8931b2e3d75a39528852bd0cb671cd16eb247f..24c3691e5a4bb932fe473558da35b269a93d4242 100644 --- a/DG_Approximation.py +++ b/DG_Approximation.py @@ -7,13 +7,15 @@ TODO: Contemplate saving 5-CV split and evaluating models separately TODO: Contemplate separating cell average and reconstruction calculations completely TODO: Contemplate removing Methods section from class docstring +TODO: Contemplate moving stiffness and boundary matrix from update scheme + to basis Urgent: TODO: Investigate protected vs. public variables -> Done TODO: Investigate decorators and properties -> Done TODO: Investigate how to make variables read-only -> Done TODO: Replace getter with property attributes for basis -> Done -TODO: Contain polynomial degree in basis +TODO: Contain polynomial degree in basis -> Done TODO: Improve docstring for Basis class TODO: Enforce abstract classes/methods (abc.ABC, abc.abstractmethod) TODO: Introduce Mesh class @@ -154,7 +156,6 @@ class DGScheme: """ # Unpack keyword arguments self._wave_speed = kwargs.pop('wave_speed', 1) - self._polynomial_degree = kwargs.pop('polynomial_degree', 2) self._cfl_number = kwargs.pop('cfl_number', 0.2) self._num_grid_cells = kwargs.pop('num_grid_cells', 64) self._final_time = kwargs.pop('final_time', 1) @@ -173,6 +174,9 @@ class DGScheme: self._quadrature_config = kwargs.pop('quadrature_config', {}) self._update_scheme = kwargs.pop('update_scheme', 'SSPRK3') + polynomial_degree = kwargs.pop('polynomial_degree', 2) + self._basis = OrthonormalLegendre(polynomial_degree) + # Throw an error if there are extra keyword arguments if len(kwargs) > 0: extra = ', '.join('"%s"' % k for k in list(kwargs.keys())) @@ -206,12 +210,11 @@ class DGScheme: self._detector = getattr(Troubled_Cell_Detector, self._detector)( config=self._detector_config, mesh=self._mesh, wave_speed=self._wave_speed, num_grid_cells=self._num_grid_cells, - polynomial_degree=self._polynomial_degree, final_time=self._final_time, left_bound=self._left_bound, right_bound=self._right_bound, basis=self._basis, init_cond=self._init_cond, quadrature=self._quadrature) self._update_scheme = getattr(Update_Scheme, self._update_scheme)( - polynomial_degree=self._polynomial_degree, + polynomial_degree=polynomial_degree, num_grid_cells=self._num_grid_cells, detector=self._detector, limiter=self._limiter) @@ -233,8 +236,7 @@ class DGScheme: projection = do_initial_projection( initial_condition=self._init_cond, basis=self._basis, quadrature=self._quadrature, num_grid_cells=self._num_grid_cells, - left_bound=self._left_bound, right_bound=self._right_bound, - polynomial_degree=self._polynomial_degree) + left_bound=self._left_bound, right_bound=self._right_bound) time_step = abs(self._cfl_number * self._cell_len / self._wave_speed) @@ -281,7 +283,6 @@ class DGScheme: # Set additional necessary instance variables self._interval_len = self._right_bound-self._left_bound self._cell_len = self._interval_len / self._num_grid_cells - self._basis = OrthonormalLegendre(self._polynomial_degree) # Set additional necessary config parameters self._limiter_config['cell_len'] = self._cell_len @@ -323,7 +324,7 @@ class DGScheme: initial_condition=initial_condition, basis=self._basis, quadrature=self._quadrature, num_grid_cells=self._num_grid_cells, left_bound=self._left_bound, right_bound=self._right_bound, - polynomial_degree=self._polynomial_degree, adjustment=adjustment) + adjustment=adjustment) return self._basis.calculate_cell_average( projection=projection[:, 1:-1], stencil_length=stencil_length, @@ -332,7 +333,7 @@ class DGScheme: def do_initial_projection(initial_condition, basis, quadrature, num_grid_cells, left_bound, right_bound, - polynomial_degree, adjustment=0): + adjustment=0): """Calculates initial projection. Calculates a projection at time step 0 and adds ghost cells on both @@ -352,8 +353,6 @@ def do_initial_projection(initial_condition, basis, quadrature, Left boundary of interval. right_bound : float Right boundary of interval. - polynomial_degree : int - Polynomial degree. adjustment: float, optional Extent of adjustment of each evaluation point in x-direction. Default: 0. @@ -362,7 +361,7 @@ def do_initial_projection(initial_condition, basis, quadrature, ------- ndarray Matrix containing projection of size (N+2, p+1) with N being the - number of grid cells and p being the polynomial degree. + number of grid cells and p being the polynomial degree of the basis. """ # Initialize matrix and set first entry to accommodate for ghost cell @@ -374,7 +373,7 @@ def do_initial_projection(initial_condition, basis, quadrature, new_row = [] eval_point = left_bound + (cell+0.5)*cell_len - for degree in range(polynomial_degree + 1): + for degree in range(basis.polynomial_degree + 1): new_row.append(np.float64(sum(initial_condition.calculate( eval_point + cell_len/2 * quadrature.get_eval_points()[point] - adjustment) diff --git a/Plotting.py b/Plotting.py index cfff3bace7357637dea9057be64bd765efe02b6b..bd104cdd9564aba07d7e6cad0b3f85ea6d9968cf 100644 --- a/Plotting.py +++ b/Plotting.py @@ -124,10 +124,9 @@ def plot_shock_tube(num_grid_cells: int, troubled_cell_history: list, plt.title('Shock Tubes') -def plot_details(fine_projection: ndarray, fine_mesh: ndarray, - coarse_projection: ndarray, basis: ndarray, wavelet: ndarray, - multiwavelet_coeffs: ndarray, num_coarse_grid_cells: int, - polynomial_degree: int) -> None: +def plot_details(fine_projection: ndarray, fine_mesh: ndarray, basis: Basis, + coarse_projection: ndarray, multiwavelet_coeffs: ndarray, + num_coarse_grid_cells: int) -> None: """Plots details of projection to coarser mesh. Parameters @@ -136,35 +135,32 @@ def plot_details(fine_projection: ndarray, fine_mesh: ndarray, Matrix of projection for each polynomial degree. fine_mesh : ndarray List of evaluation points for fine mesh. - basis : ndarray - Basis vector for calculation. - wavelet : ndarray - Wavelet vector for calculation. + basis: Basis object + Basis used for calculation. multiwavelet_coeffs : ndarray Matrix of multiwavelet coefficients. num_coarse_grid_cells : int Number of cells in the coarse mesh (half the cells of the fine mesh). Usually exponential of 2. - polynomial_degree : int - Polynomial degree. """ averaged_projection = [[coarse_projection[degree][cell] - * basis[degree].subs(x, value) + * basis.basis[degree].subs(x, value) for cell in range(num_coarse_grid_cells) for value in [-0.5, 0.5]] - for degree in range(polynomial_degree + 1)] + for degree in range(basis.polynomial_degree + 1)] wavelet_projection = [[multiwavelet_coeffs[degree][cell] - * wavelet[degree].subs(z, 0.5) * value + * basis.wavelet[degree].subs(z, 0.5) * value for cell in range(num_coarse_grid_cells) - for value in [(-1) ** (polynomial_degree + for value in [(-1) ** (basis.polynomial_degree + degree + 1), 1]] - for degree in range(polynomial_degree + 1)] + for degree in range(basis.polynomial_degree + 1)] projected_coarse = np.sum(averaged_projection, axis=0) - projected_fine = np.sum([fine_projection[degree] * basis[degree].subs(x, 0) - for degree in range(polynomial_degree + 1)], + projected_fine = np.sum([fine_projection[degree] + * basis.basis[degree].subs(x, 0) + for degree in range(basis.polynomial_degree + 1)], axis=0) projected_wavelet_coeffs = np.sum(wavelet_projection, axis=0) @@ -352,6 +348,7 @@ def plot_approximation_results(data_file: str, directory: str, plot_name: str, # Decode all ndarrays by converting lists approx_stats = {key: decode_ndarray(approx_stats[key]) for key in approx_stats.keys()} + approx_stats.pop('polynomial_degree') # Plot exact/approximate results, errors, shock tubes, # and any detector-dependant plots @@ -375,7 +372,7 @@ def plot_approximation_results(data_file: str, directory: str, plot_name: str, def plot_results(projection: ndarray, troubled_cell_history: list, time_history: list, mesh: ndarray, num_grid_cells: int, - polynomial_degree: int, wave_speed: float, final_time: float, + wave_speed: float, final_time: float, left_bound: float, right_bound: float, basis: Basis, quadrature: Quadrature, init_cond: InitialCondition, colors: dict = None, coarse_projection: ndarray = None, @@ -400,8 +397,6 @@ def plot_results(projection: ndarray, troubled_cell_history: list, List of mesh valuation points. num_grid_cells : int Number of cells in the mesh. Usually exponential of 2. - polynomial_degree : int - Polynomial degree. wave_speed : float Speed of wave in rightward direction. final_time : float @@ -444,7 +439,7 @@ def plot_results(projection: ndarray, troubled_cell_history: list, init_cond) approx = calculate_approximate_solution( projection[:, 1:-1], quadrature.get_eval_points(), - polynomial_degree, basis.basis) + basis.polynomial_degree, basis.basis) # Plot multiwavelet solution (fine and coarse grid) if coarse_projection is not None: @@ -460,17 +455,15 @@ def plot_results(projection: ndarray, troubled_cell_history: list, init_cond) coarse_approx = calculate_approximate_solution( coarse_projection, quadrature.get_eval_points(), - polynomial_degree, basis.basis) + basis.polynomial_degree, basis.basis) plot_solution_and_approx( coarse_grid, coarse_exact, coarse_approx, colors['coarse_exact'], colors['coarse_approx']) # Plot multiwavelet details num_coarse_grid_cells = num_grid_cells//2 - plot_details(projection[:, 1:-1], mesh[2:-2], coarse_projection, - basis.basis, basis.wavelet, multiwavelet_coeffs, - num_coarse_grid_cells, - polynomial_degree) + plot_details(projection[:, 1:-1], mesh[2:-2], basis, coarse_projection, + multiwavelet_coeffs, num_coarse_grid_cells) plot_solution_and_approx(grid, exact, approx, colors['fine_exact'], @@ -491,7 +484,7 @@ def plot_results(projection: ndarray, troubled_cell_history: list, plot_semilog_error(grid, pointwise_error) plot_error(grid, exact, approx) - print('p =', polynomial_degree) + print('p =', basis.polynomial_degree) print('N =', num_grid_cells) print('maximum error =', max_error) diff --git a/Troubled_Cell_Detector.py b/Troubled_Cell_Detector.py index 6db101be3b8ef1d7f5e7f9d10c96b3a3767bd4a6..f73caefbb94d8b8992ae1a01c9432f7c280f265b 100644 --- a/Troubled_Cell_Detector.py +++ b/Troubled_Cell_Detector.py @@ -37,7 +37,7 @@ class TroubledCellDetector: """ def __init__(self, config, init_cond, quadrature, basis, mesh, - wave_speed=1, polynomial_degree=2, num_grid_cells=64, + wave_speed=1, num_grid_cells=64, final_time=1, left_bound=-1, right_bound=1): """Initializes TroubledCellDetector. @@ -55,8 +55,6 @@ class TroubledCellDetector: List of mesh valuation points. wave_speed : float, optional Speed of wave in rightward direction. Default: 1. - polynomial_degree : int, optional - Polynomial degree. Default: 2. num_grid_cells : int, optional Number of cells in the mesh. Usually exponential of 2. Default: 64. final_time : float, optional @@ -69,7 +67,6 @@ class TroubledCellDetector: """ self._mesh = mesh self._wave_speed = wave_speed - self._polynomial_degree = polynomial_degree self._num_grid_cells = num_grid_cells self._final_time = final_time self._left_bound = left_bound @@ -113,7 +110,7 @@ class TroubledCellDetector: 'num_grid_cells': self._num_grid_cells, 'mesh': self._mesh, 'final_time': self._final_time, 'left_bound': self._left_bound, 'right_bound': self._right_bound, - 'polynomial_degree': self._polynomial_degree + 'polynomial_degree': self._basis.polynomial_degree } @@ -513,7 +510,8 @@ class Theoretical(WaveletDetector): """ max_value = max(abs(multiwavelet_coeffs[degree][cell]) - for degree in range(self._polynomial_degree+1))/max_avg + for degree in range( + self._basis.polynomial_degree+1))/max_avg eps = self._cutoff_factor\ / (self._cell_len*self._num_coarse_grid_cells*2)