From 2c7ae512a6650821109d1feb63f036f5f0ed62bd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?K=C3=BChle=2C=20Laura=20Christine=20=28lakue103=29?= <laura.kuehle@uni-duesseldorf.de> Date: Thu, 24 Nov 2022 20:38:31 +0100 Subject: [PATCH] Changed code to allow compatibility with arbitrary number of ghost cells. --- scripts/approximate_solution.py | 11 +++---- scripts/tcd/ANN_Data_Generator.py | 3 +- scripts/tcd/DG_Approximation.py | 13 +++++--- scripts/tcd/Limiter.py | 33 +++++++++++++------- scripts/tcd/Mesh.py | 5 +-- scripts/tcd/Plotting.py | 8 +++-- scripts/tcd/Troubled_Cell_Detector.py | 15 ++++++--- scripts/tcd/Update_Scheme.py | 44 ++++++++++++++++++--------- 8 files changed, 85 insertions(+), 47 deletions(-) diff --git a/scripts/approximate_solution.py b/scripts/approximate_solution.py index 32fc785..bbda5a4 100644 --- a/scripts/approximate_solution.py +++ b/scripts/approximate_solution.py @@ -35,7 +35,7 @@ def main() -> None: mesh = Mesh(num_cells=params.pop('num_mesh_cells', 64), left_bound=params.pop('left_bound', -1), right_bound=params.pop('right_bound', 1), - num_ghost_cells=2) + num_ghost_cells=1) # Initialize basis basis = OrthonormalLegendre(params.pop('polynomial_degree', 2)) @@ -43,10 +43,10 @@ def main() -> None: # Initialize limiter after checking its legitimacy limiter_name = params.pop('limiter', 'ModifiedMinMod') limiter_config = params.pop('limiter_config', {}) - limiter_config['cell_len'] = mesh.cell_len if not hasattr(Limiter, limiter_name): raise ValueError('Invalid limiter: "%s"' % limiter_name) - limiter = getattr(Limiter, limiter_name)(config=limiter_config) + limiter = getattr(Limiter, limiter_name)(mesh=mesh, + config=limiter_config) # Initialize troubled cell detector after checking its legitimacy detector_name = params.pop('detector') @@ -63,9 +63,8 @@ def main() -> None: if not hasattr(Update_Scheme, scheme_name): raise ValueError('Invalid update scheme: "%s"' % scheme_name) update_scheme = getattr(Update_Scheme, scheme_name)( - polynomial_degree=basis.polynomial_degree, - num_cells=mesh.num_cells, detector=detector, - limiter=limiter, wave_speed=wave_speed) + polynomial_degree=basis.polynomial_degree, mesh=mesh, + detector=detector, limiter=limiter, wave_speed=wave_speed) # Initialize quadrature after checking its legitimacy quadrature_name = params.pop('quadrature', 'Gauss') diff --git a/scripts/tcd/ANN_Data_Generator.py b/scripts/tcd/ANN_Data_Generator.py index 184a93d..881f4d2 100644 --- a/scripts/tcd/ANN_Data_Generator.py +++ b/scripts/tcd/ANN_Data_Generator.py @@ -226,8 +226,7 @@ class TrainingDataGenerator: x_shift=shift) input_data[i] = self._basis_list[ polynomial_degree].calculate_cell_average( - projection=projection[:, 1:-1], - stencil_len=stencil_len, + projection=projection, stencil_len=stencil_len, add_reconstructions=add_reconstructions) count += 1 diff --git a/scripts/tcd/DG_Approximation.py b/scripts/tcd/DG_Approximation.py index a816f08..7df9b44 100644 --- a/scripts/tcd/DG_Approximation.py +++ b/scripts/tcd/DG_Approximation.py @@ -182,7 +182,8 @@ def do_initial_projection(init_cond, mesh, basis, quadrature, """ # Initialize output matrix - output_matrix = np.zeros((basis.polynomial_degree+1, mesh.num_cells+2)) + output_matrix = np.zeros((basis.polynomial_degree+1, + mesh.num_cells+2*mesh.num_ghost_cells)) # Calculate basis based on quadrature basis_matrix = np.array([np.vectorize(basis.basis[degree].subs)( @@ -196,10 +197,14 @@ def do_initial_projection(init_cond, mesh, basis, quadrature, x=points, mesh=mesh) # Set output matrix for regular cells - output_matrix[:, 1:-1] = (basis_matrix * quadrature.weights) @ init_matrix + output_matrix[:, mesh.num_ghost_cells: + output_matrix.shape[-1]-mesh.num_ghost_cells] = \ + (basis_matrix * quadrature.weights) @ init_matrix # Set ghost cells - output_matrix[:, 0] = output_matrix[:, -2] - output_matrix[:, -1] = output_matrix[:, 1] + output_matrix[:, :mesh.num_ghost_cells] = \ + output_matrix[:, -2*mesh.num_ghost_cells:-mesh.num_ghost_cells] + output_matrix[:, output_matrix.shape[-1]-mesh.num_ghost_cells:] = \ + output_matrix[:, mesh.num_ghost_cells:2*mesh.num_ghost_cells] return output_matrix diff --git a/scripts/tcd/Limiter.py b/scripts/tcd/Limiter.py index 258cee1..3042c55 100644 --- a/scripts/tcd/Limiter.py +++ b/scripts/tcd/Limiter.py @@ -18,15 +18,19 @@ class Limiter(ABC): Applies limiting to cells. """ - def __init__(self, config): + def __init__(self, mesh, config): """Initializes Quadrature. Parameters ---------- + mesh : Mesh + Mesh for calculation. config : dict Additional parameters for limiter. """ + self._mesh = mesh + self._reset(config) @abstractmethod @@ -144,9 +148,9 @@ class MinMod(Limiter): modification_mask = self._determine_mask(projection, cell_slopes) cells = np.array(cells) mask = np.zeros_like(new_projection, dtype=bool) - mask[self._erase_degree+1:, cells+1] = np.tile( - np.logical_not(modification_mask)[cells], - (len(projection)-self._erase_degree-1, 1)) + mask[self._erase_degree+1:, cells+self._mesh.num_ghost_cells] = \ + np.tile(np.logical_not(modification_mask)[cells], + (len(projection)-self._erase_degree-1, 1)) # Limit troubled cells for higher degrees new_projection[mask] = 0 @@ -169,8 +173,17 @@ class MinMod(Limiter): Mask whether cells should be adjusted. """ - forward_slopes = (projection[0, 2:]-projection[0, 1:-1]) * (0.5**0.5) - backward_slopes = (projection[0, 1:-1]-projection[0, :-2]) * (0.5**0.5) + forward_slopes = (projection[0, self._mesh.num_ghost_cells+1: + projection.shape[-1] + - self._mesh.num_ghost_cells+1] - + projection[0, self._mesh.num_ghost_cells: + -self._mesh.num_ghost_cells]) * \ + (0.5**0.5) + backward_slopes = (projection[0, self._mesh.num_ghost_cells: + -self._mesh.num_ghost_cells] - + projection[0, self._mesh.num_ghost_cells-1: + - self._mesh.num_ghost_cells-1]) * \ + (0.5**0.5) pos_mask = np.logical_and(slopes >= 0, np.logical_and(forward_slopes >= 0, backward_slopes >= 0)) @@ -181,8 +194,7 @@ class MinMod(Limiter): return slope_mask - @staticmethod - def _set_cell_slope(projection): + def _set_cell_slope(self, projection): """Calculates the slope of the cell. Parameters @@ -199,7 +211,7 @@ class MinMod(Limiter): root_vector = np.array([np.sqrt(degree+0.5) for degree in range(len(projection))]) slope = root_vector[1:] @ projection[1:] - return slope[1:-1] + return slope[self._mesh.num_ghost_cells:-self._mesh.num_ghost_cells] class ModifiedMinMod(MinMod): @@ -236,9 +248,8 @@ class ModifiedMinMod(MinMod): super()._reset(config) # Unpack necessary configurations - cell_len = config.pop('cell_len') mod_factor = config.pop('mod_factor', 0) - self._threshold = mod_factor * cell_len**2 + self._threshold = mod_factor * self._mesh.cell_len**2 def get_name(self): """Returns string of class name concatenated with the erase-degree.""" diff --git a/scripts/tcd/Mesh.py b/scripts/tcd/Mesh.py index 36abe12..4f0a526 100644 --- a/scripts/tcd/Mesh.py +++ b/scripts/tcd/Mesh.py @@ -121,7 +121,8 @@ class Mesh: @property def non_ghost_cells(self) -> ndarray: """Return the cell centers of the mesh (excluding ghost cells).""" - return self.cells[self._num_ghost_cells:-self._num_ghost_cells] + return self.cells[self._num_ghost_cells: + len(self.cells)-self._num_ghost_cells] def create_data_dict(self) -> dict: """Return dictionary with data necessary to construct mesh.""" @@ -155,5 +156,5 @@ class Mesh: # Return new mesh instance return Mesh(left_bound=point - stencil_len/2 * mesh_spacing, right_bound=point + stencil_len/2 * mesh_spacing, - num_cells=stencil_len, num_ghost_cells=2, + num_cells=stencil_len, num_ghost_cells=0, training_data_mode=True) diff --git a/scripts/tcd/Plotting.py b/scripts/tcd/Plotting.py index 14ffcf4..42f8fdc 100644 --- a/scripts/tcd/Plotting.py +++ b/scripts/tcd/Plotting.py @@ -419,14 +419,16 @@ def plot_results(projection: ndarray, troubled_cell_history: list, # Determine exact and approximate solution grid, exact = calculate_exact_solution( mesh, wave_speed, final_time, quadrature, init_cond) + projection = projection[:, mesh.num_ghost_cells: + -mesh.num_ghost_cells] approx = calculate_approximate_solution( - projection[:, 1:-1], quadrature.nodes, + projection, quadrature.nodes, basis.polynomial_degree, basis.basis) # Plot multiwavelet solution (fine and coarse mesh) if coarse_projection is not None: coarse_mesh = Mesh(num_cells=mesh.num_cells//2, - num_ghost_cells=1, left_bound=mesh.bounds[0], + num_ghost_cells=0, left_bound=mesh.bounds[0], right_bound=mesh.bounds[1]) # Plot exact and approximate solutions for coarse mesh @@ -440,7 +442,7 @@ def plot_results(projection: ndarray, troubled_cell_history: list, colors['coarse_approx']) # Plot multiwavelet details - plot_details(projection[:, 1:-1], mesh, basis, coarse_projection, + plot_details(projection, mesh, basis, coarse_projection, multiwavelet_coeffs) plot_solution_and_approx(grid, exact, approx, diff --git a/scripts/tcd/Troubled_Cell_Detector.py b/scripts/tcd/Troubled_Cell_Detector.py index 175fd9c..805fdd8 100644 --- a/scripts/tcd/Troubled_Cell_Detector.py +++ b/scripts/tcd/Troubled_Cell_Detector.py @@ -175,7 +175,8 @@ class ArtificialNeuralNetwork(TroubledCellDetector): """ # Reset ghost cells to adjust for stencil length num_ghost_cells = self._stencil_len//2 - projection = projection[:, 1:-1] + projection = projection[:, self._mesh.num_ghost_cells: + -self._mesh.num_ghost_cells] projection = np.concatenate((projection[:, -num_ghost_cells:], projection, projection[:, :num_ghost_cells]), axis=1) @@ -238,8 +239,13 @@ class WaveletDetector(TroubledCellDetector): Matrix of multiwavelet coefficients. """ - return 0.5 * (self._wavelet_projection_left.T @ projection[:, 1:-1] + - self._wavelet_projection_right.T @ projection[:, 2:]) + return 0.5 * (self._wavelet_projection_left.T @ + projection[:, self._mesh.num_ghost_cells: + -self._mesh.num_ghost_cells] + + self._wavelet_projection_right.T @ + projection[:, self._mesh.num_ghost_cells+1: + projection.shape[-1] + - self._mesh.num_ghost_cells+1]) def _select_degree(self, wavelet_matrix): """Select degree of wavelet coefficients for troubled cell detection. @@ -285,7 +291,8 @@ class WaveletDetector(TroubledCellDetector): = self._basis.basis_projection # Remove ghost cells - projection = projection[:, 1:-1] + projection = projection[:, self._mesh.num_ghost_cells: + -self._mesh.num_ghost_cells] # Calculate projection on coarse mesh return 0.5 * (basis_projection_left.T @ projection[:, ::2] + diff --git a/scripts/tcd/Update_Scheme.py b/scripts/tcd/Update_Scheme.py index b418732..19cd201 100644 --- a/scripts/tcd/Update_Scheme.py +++ b/scripts/tcd/Update_Scheme.py @@ -26,7 +26,7 @@ class UpdateScheme(ABC): Performs time step. """ - def __init__(self, polynomial_degree, num_cells, detector, limiter, + def __init__(self, polynomial_degree, mesh, detector, limiter, wave_speed): """Initializes UpdateScheme. @@ -34,8 +34,8 @@ class UpdateScheme(ABC): ---------- polynomial_degree : int Polynomial degree. - num_cells : int - Number of cells in the mesh. Usually exponential of 2. + mesh : Mesh + Mesh for calculation. detector : TroubledCellDetector object Troubled cell detector for evaluation. limiter : Limiter object @@ -46,7 +46,7 @@ class UpdateScheme(ABC): """ # Unpack positional arguments self._polynomial_degree = polynomial_degree - self._num_cells = num_cells + self._mesh = mesh self._detector = detector self._limiter = limiter self._wave_speed = wave_speed @@ -168,9 +168,12 @@ class UpdateScheme(ABC): degree. """ - current_projection[:, 0] = current_projection[:, self._num_cells] - current_projection[:, self._num_cells+1] \ - = current_projection[:, 1] + current_projection[:, :self._mesh.num_ghost_cells] = \ + current_projection[ + :, -2*self._mesh.num_ghost_cells:-self._mesh.num_ghost_cells] + current_projection[:, -self._mesh.num_ghost_cells:] = \ + current_projection[ + :, self._mesh.num_ghost_cells:2*self._mesh.num_ghost_cells] return current_projection @@ -311,16 +314,27 @@ class SSPRK3(UpdateScheme): """ right_hand_side = np.zeros_like(current_projection) if self._wave_speed > 0: - right_hand_side[:, 1:-1] = 2 * ( - self._boundary_matrix @ current_projection[:, :-2] + - self._stiffness_matrix @ current_projection[:, 1:-1]) + right_hand_side[:, self._mesh.num_ghost_cells: + -self._mesh.num_ghost_cells] = \ + 2 * (self._boundary_matrix @ + current_projection[:, self._mesh.num_ghost_cells-1: + -self._mesh.num_ghost_cells-1] + + self._stiffness_matrix @ + current_projection[:, self._mesh.num_ghost_cells: + -self._mesh.num_ghost_cells]) else: - right_hand_side[:, 1:-1] = 2 * ( - self._boundary_matrix @ current_projection[:, 2:] + - self._stiffness_matrix @ current_projection[:, 1:-1]) + right_hand_side[:, self._mesh.num_ghost_cells: + -self._mesh.num_ghost_cells] = \ + 2 * (self._boundary_matrix @ + current_projection[:, self._mesh.num_ghost_cells+1:] + + self._stiffness_matrix @ + current_projection[:, self._mesh.num_ghost_cells: + -self._mesh.num_ghost_cells]) # Set ghost cells to respective value - right_hand_side[:, 0] = right_hand_side[:, -2] - right_hand_side[:, -1] = right_hand_side[:, 1] + right_hand_side[:, :self._mesh.num_ghost_cells] = right_hand_side[ + :, -2*self._mesh.num_ghost_cells:-self._mesh.num_ghost_cells] + right_hand_side[:, -self._mesh.num_ghost_cells:] = right_hand_side[ + :, self._mesh.num_ghost_cells:2*self._mesh.num_ghost_cells] return right_hand_side -- GitLab