diff --git a/Snakefile b/Snakefile index 615f4a759d04993926f8088426838b039a541ee7..6cd65c6fb45372d8c5d22eed0e1f118ad47dd44c 100644 --- a/Snakefile +++ b/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 diff --git a/scripts/tcd/Boundary_Condition.py b/scripts/tcd/Boundary_Condition.py new file mode 100644 index 0000000000000000000000000000000000000000..3d5ab91bcaaaa4c1960b3118e99cd0fc10b5cfb6 --- /dev/null +++ b/scripts/tcd/Boundary_Condition.py @@ -0,0 +1,76 @@ +# -*- 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 diff --git a/scripts/tcd/Equation.py b/scripts/tcd/Equation.py index a4513ed5898fbc82c43415d50da1655bcb4b2d6c..62c9ba32463770813780f1db504264501b2e7b3e 100644 --- a/scripts/tcd/Equation.py +++ b/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 diff --git a/scripts/tcd/Limiter.py b/scripts/tcd/Limiter.py index 3042c55afbfe3a128d2debd8f8e14950376994ac..7fcec4c1e1d1ba2e794c96a438e594c1ac521501 100644 --- a/scripts/tcd/Limiter.py +++ b/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. diff --git a/scripts/tcd/Mesh.py b/scripts/tcd/Mesh.py index 5d5d1d61183ad1fb94f8fdbf2a80ed1cf319093c..649fd55799dbc9287b1fb3fa7a5145a16b7761d9 100644 --- a/scripts/tcd/Mesh.py +++ b/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: diff --git a/scripts/tcd/Update_Scheme.py b/scripts/tcd/Update_Scheme.py index e7c9122b13803a0d41fecd91bfac7177ebd21abe..ee83bac7fe0893349717e152caefddc4514137c9 100644 --- a/scripts/tcd/Update_Scheme.py +++ b/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 diff --git a/scripts/tcd/projection_utils.py b/scripts/tcd/projection_utils.py index 17e0b3e1dcde2d33ac1bf0d86868fabc11d361d0..202f50e199d34a6252568272d44d07670c196999 100644 --- a/scripts/tcd/projection_utils.py +++ b/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