Changed code to allow compatibility with arbitrary number of ghost cells.

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')

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

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

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,8 +148,8 @@ 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],
+        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
@@ -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."""

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)

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,

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] +

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