Enforced periodic boundary condition for projection with decorator.

Snakefile

@@ -24,7 +24,7 @@ TODO: Discuss how wavelet details should be plotted
 
 Urgent:
 TODO: Refactor Boxplot class -> Done
-TODO: Enforce boundary condition with decorator
+TODO: Enforce periodic boundary condition for projection with decorator -> Done
 TODO: Enforce Boxplot bounds with decorator
 TODO: Enforce Boxplot folds with decorator
 TODO: Enforce boundary for initial condition in exact solution only

scripts/tcd/Boundary_Condition.py

+# -*- coding: utf-8 -*-
+"""
+@author: Laura C. Kühle
+
+"""
+# from typing import Tuple
+# from abc import ABC, abstractmethod
+# import numpy as np
+from numpy import ndarray
+
+# from .Basis_Function import Basis
+# from .Initial_Condition import InitialCondition
+# from .Mesh import Mesh
+# from .projection_utils import do_initial_projection
+# from .Quadrature import Quadrature
+
+
+def periodic_boundary(projection: ndarray, num_ghost_cells: int) -> ndarray:
+    """Enforce boundary condition.
+
+    Adjust ghost cells to ensure periodic boundary condition.
+
+    Parameters
+    ----------
+    projection : ndarray
+        Matrix of projection for each polynomial degree.
+    num_ghost_cells : int
+        Number of ghost cells to be adjusted.
+
+    Returns
+    -------
+    projection : ndarray
+        Matrix of projection for each polynomial degree.
+
+    """
+    projection[:, :num_ghost_cells] = \
+        projection[:, -2 * num_ghost_cells:-num_ghost_cells]
+    projection[:, -num_ghost_cells:] = \
+        projection[:, num_ghost_cells:2 * num_ghost_cells]
+    return projection
+
+
+def dirichlet_boundary(projection: ndarray, num_ghost_cells: int) -> ndarray:
+    """Enforce boundary condition.
+
+    Adjust ghost cells to ensure Dirichlet boundary condition.
+
+    Parameters
+    ----------
+    projection : ndarray
+        Matrix of projection for each polynomial degree.
+    num_ghost_cells : int
+        Number of ghost cells to be adjusted.
+
+    Returns
+    -------
+    projection : ndarray
+        Matrix of projection for each polynomial degree.
+
+    """
+    pass
+
+
+def enforce_boundary(num_ghost_cells=None):
+    def _enforce_boundary(func):
+        def boundary(self, *args, **kwargs):
+            projection = func(self, *args, **kwargs)
+            if self._mesh.boundary == 'periodic':
+                return periodic_boundary(
+                    projection=projection, num_ghost_cells=num_ghost_cells
+                    if num_ghost_cells is not None
+                    else self._mesh.num_ghost_cells)
+            else:
+                raise Exception('Not implemented!')
+        return boundary
+    return _enforce_boundary

scripts/tcd/Equation.py

@@ -9,6 +9,7 @@ import numpy as np
 from numpy import ndarray
 
 from .Basis_Function import Basis
+from .Boundary_Condition import enforce_boundary
 from .Initial_Condition import InitialCondition
 from .Mesh import Mesh
 from .projection_utils import do_initial_projection
@@ -26,8 +27,6 @@ class Equation(ABC):
         Update time step.
     solve_exactly(mesh)
         Calculate exact solution.
-    enforce_boundary_condition(projection)
-        Enforce boundary condition.
     update_right_hand_side(projection)
         Update right-hand side.
     """
@@ -93,6 +92,7 @@ class Equation(ABC):
         return abs(self._cfl_number * self._mesh.cell_len / self._wave_speed)
 
     @abstractmethod
+    @enforce_boundary()
     def _initialize_projection(self) -> ndarray:
         """Initialize projection."""
         pass
@@ -138,25 +138,7 @@ class Equation(ABC):
         pass
 
     @abstractmethod
-    def enforce_boundary_condition(self, projection: ndarray) -> ndarray:
-        """Enforce boundary condition.
-
-        Adjust ghost cells to ensure periodic boundary condition.
-
-        Parameters
-        ----------
-        projection : ndarray
-            Matrix of projection for each polynomial degree.
-
-        Returns
-        -------
-        projection : ndarray
-            Matrix of projection for each polynomial degree.
-
-        """
-        pass
-
-    @abstractmethod
+    @enforce_boundary()
     def update_right_hand_side(self, projection: ndarray) -> ndarray:
         """Update right-hand side.
 
@@ -188,8 +170,6 @@ class LinearAdvection(Equation):
     -------
     update_right_hand_side(projection)
         Update right-hand side.
-    enforce_boundary_condition(projection)
-        Enforce boundary condition.
 
     Notes
     -----
@@ -222,6 +202,7 @@ class LinearAdvection(Equation):
             else np.fromfunction(lambda i, j: (-1.0)**j, matrix_shape)
         self._flux_matrix = matrix * degree_matrix
 
+    @enforce_boundary()
     def _initialize_projection(self) -> ndarray:
         """Initialize projection."""
         return do_initial_projection(init_cond=self._init_cond,
@@ -283,28 +264,7 @@ class LinearAdvection(Equation):
 
         return grid, exact
 
-    def enforce_boundary_condition(self, projection: ndarray) -> ndarray:
-        """Enforce boundary condition.
-
-        Adjust ghost cells to ensure periodic boundary condition.
-
-        Parameters
-        ----------
-        projection : ndarray
-            Matrix of projection for each polynomial degree.
-
-        Returns
-        -------
-        projection : ndarray
-            Matrix of projection for each polynomial degree.
-
-        """
-        projection[:, :self._mesh.num_ghost_cells] = projection[
-            :, -2*self._mesh.num_ghost_cells:-self._mesh.num_ghost_cells]
-        projection[:, -self._mesh.num_ghost_cells:] = projection[
-            :, self._mesh.num_ghost_cells:2*self._mesh.num_ghost_cells]
-        return projection
-
+    @enforce_boundary()
     def update_right_hand_side(self, projection: ndarray) -> ndarray:
         """Update right-hand side.
 
@@ -338,10 +298,4 @@ class LinearAdvection(Equation):
                      projection[:, self._mesh.num_ghost_cells:
                                 -self._mesh.num_ghost_cells])
 
-        # Set ghost cells to respective value
-        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

scripts/tcd/Limiter.py

@@ -6,6 +6,8 @@
 from abc import ABC, abstractmethod
 import numpy as np
 
+from .Boundary_Condition import enforce_boundary
+
 
 class Limiter(ABC):
     """Abstract class for limiting function.
@@ -50,6 +52,7 @@ class Limiter(ABC):
         return self.__class__.__name__
 
     @abstractmethod
+    @enforce_boundary()
     def apply(self, projection, cells):
         """Applies limiting to cells.
 
@@ -80,6 +83,7 @@ class NoLimiter(Limiter):
         Applies no limiting to cells.
 
     """
+    @enforce_boundary()
     def apply(self, projection, cells):
         """Returns projection without limiting."""
         return projection
@@ -121,6 +125,7 @@ class MinMod(Limiter):
         """Returns string of class name concatenated with the erase-degree."""
         return self.__class__.__name__ + str(self._erase_degree)
 
+    @enforce_boundary()
     def apply(self, projection, cells):
         """Applies limiting to cells.
 

scripts/tcd/Mesh.py

@@ -77,6 +77,7 @@ class Mesh:
                                  'an exponential of 2')
         self._left_bound = left_bound
         self._right_bound = right_bound
+        self.boundary = 'periodic'
 
     @property
     def mode(self) -> str:

scripts/tcd/Update_Scheme.py

@@ -150,23 +150,17 @@ class SSPRK3(UpdateScheme):
         current_projection = self._apply_first_step(original_projection,
                                                     cfl_number)
         current_projection, __ = self._apply_limiter(current_projection)
-        current_projection = self._equation.enforce_boundary_condition(
-            current_projection)
 
         current_projection = self._apply_second_step(original_projection,
                                                      current_projection,
                                                      cfl_number)
         current_projection, __ = self._apply_limiter(current_projection)
-        current_projection = self._equation.enforce_boundary_condition(
-            current_projection)
 
         current_projection = self._apply_third_step(original_projection,
                                                     current_projection,
                                                     cfl_number)
         current_projection, troubled_cells = self._apply_limiter(
             current_projection)
-        current_projection = self._equation.enforce_boundary_condition(
-            current_projection)
 
         return current_projection, troubled_cells
 

scripts/tcd/projection_utils.py

@@ -140,10 +140,4 @@ def do_initial_projection(init_cond, mesh, basis, quadrature,
                   output_matrix.shape[-1]-mesh.num_ghost_cells] = \
         (basis_matrix * quadrature.weights) @ init_matrix
 
-    # Set ghost cells
-    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