diff --git a/Snakefile b/Snakefile
index 03a1175b4c6bacaab0fee29652dffc9b78d8bc21..6616d907e57a90435b3435f612ed82c5516ad1fd 100644
--- a/Snakefile
+++ b/Snakefile
@@ -24,7 +24,8 @@ TODO: Discuss how wavelet details should be plotted
TODO: Ask whether we are dealing with inviscid Burgers only
TODO: Ask whether it is helpful to enforce boundary at every SSPRK3 step
TODO: Ask why we limit the initial time step
-TODO: ASk whether cfl number should be absolute value
+TODO: Ask whether cfl number should be absolute value
+TODO: Ask why cells 1 and -2 are adjusted for Dirichlet boundary
Urgent:
TODO: Add Burgers class completely
@@ -33,7 +34,9 @@ TODO: Add derivative basis to Basis class -> Done
TODO: Fix wrong type for number of quadrature nodes -> Done
TODO: Add 'solve_exactly()' to Burgers
TODO: Add 'update_right_hand_side()' to Burgers -> Done
-TODO: Add Dirichlet boundary
+TODO: Add Dirichlet boundary -> Done
+TODO: Add functions to create each object from dict
+TODO: Initialize boundary config properly
TODO: Test Burgers
TODO: Rework Burgers for speed-up and better structure
TODO: Add initialization to Burgers
@@ -63,7 +66,6 @@ TODO: Add tests for all functions
Not feasible yet or doc-related:
TODO: Extract objects from UpdateScheme
TODO: Clean up result plotting (remove object properties of Equation)
-TODO: Add functions to create each object from dict
TODO: Move plot_approximation_results() into plotting script
TODO: Move plot_results() into plotting script
TODO: Move plot_evaluation_results() into plotting script
diff --git a/scripts/tcd/Basis_Function.py b/scripts/tcd/Basis_Function.py
index 69b134701c675316008201abc37712b521a57fcb..4068a8e3682c216e6f45bade8aa323e6e2fa931e 100644
--- a/scripts/tcd/Basis_Function.py
+++ b/scripts/tcd/Basis_Function.py
@@ -22,7 +22,7 @@ class Basis(ABC):
Attributes
----------
- polynomial degree : int
+ polynomial_degree : int
Polynomial degree.
basis : ndarray
Array of basis.
diff --git a/scripts/tcd/Boundary_Condition.py b/scripts/tcd/Boundary_Condition.py
index 8375491a0c4efbb3a39f52dcf43df9c0a284ac8d..1cde4ea28cb6c6f839a0564c6da446c11bdc1348 100644
--- a/scripts/tcd/Boundary_Condition.py
+++ b/scripts/tcd/Boundary_Condition.py
@@ -7,7 +7,7 @@ import numpy as np
from numpy import ndarray
-def periodic_boundary(projection: ndarray, num_ghost_cells: int) -> ndarray:
+def periodic_boundary(projection: ndarray, config: dict) -> ndarray:
"""Enforce boundary condition.
Adjust ghost cells to ensure periodic boundary condition.
@@ -16,8 +16,8 @@ def periodic_boundary(projection: ndarray, num_ghost_cells: int) -> ndarray:
----------
projection : ndarray
Matrix of projection for each polynomial degree.
- num_ghost_cells : int
- Number of ghost cells to be adjusted.
+ config : dict
+ Configuration dictionary.
Returns
-------
@@ -25,6 +25,8 @@ def periodic_boundary(projection: ndarray, num_ghost_cells: int) -> ndarray:
Matrix of projection for each polynomial degree.
"""
+ num_ghost_cells = config['num_ghost_cells']
+
projection[:, :num_ghost_cells] = \
projection[:, -2 * num_ghost_cells:-num_ghost_cells]
projection[:, -num_ghost_cells:] = \
@@ -32,7 +34,7 @@ def periodic_boundary(projection: ndarray, num_ghost_cells: int) -> ndarray:
return projection
-def dirichlet_boundary(projection: ndarray, num_ghost_cells: int) -> ndarray:
+def dirichlet_boundary(projection: ndarray, config: dict) -> ndarray:
"""Enforce boundary condition.
Adjust ghost cells to ensure Dirichlet boundary condition.
@@ -41,8 +43,8 @@ def dirichlet_boundary(projection: ndarray, num_ghost_cells: int) -> ndarray:
----------
projection : ndarray
Matrix of projection for each polynomial degree.
- num_ghost_cells : int
- Number of ghost cells to be adjusted.
+ config : dict
+ Configuration dictionary.
Returns
-------
@@ -50,18 +52,43 @@ def dirichlet_boundary(projection: ndarray, num_ghost_cells: int) -> ndarray:
Matrix of projection for each polynomial degree.
"""
- pass
+ polynomial_degree = config['polynomial_degree']
+ final_time = config['final_time']
+ num_ghost_cells = config['num_ghost_cells']
+ left_factor = config['left_factor']
+ right_factor = config['right_factor']
+
+ projection_values_left = [0 for i in range(polynomial_degree+1)]
+ projection_values_right = [0 for i in range(polynomial_degree+1)]
+ if final_time > 1:
+ projection_values_left[0] = -np.sqrt(2) / final_time
+ projection_values_right[0] = np.sqrt(2) / final_time
+ else:
+ projection_values_left[0] = np.sqrt(2) * left_factor
+ projection_values_right[0] = np.sqrt(2) * right_factor
+
+ projection[:, :num_ghost_cells+1] = np.repeat(projection_values_left)
+ projection[:, -num_ghost_cells-1:] = np.repeat(projection_values_right)
+ # projection[:, 0] = projection_values_left
+ # projection[:, 1] = projection_values_left
+ # projection[:, -2] = projection_values_right
+ # projection[:, -1] = projection_values_right
+ return projection
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':
+ if self._mesh.boundary_config['type'] == 'periodic':
+ config = self._mesh.boundary_config.copy()
+ if num_ghost_cells is not None:
+ config['num_ghost_cells'] = num_ghost_cells
return periodic_boundary(
- projection=projection, num_ghost_cells=num_ghost_cells
- if num_ghost_cells is not None
- else self._mesh.num_ghost_cells)
+ projection=projection, config=config)
+ elif self._mesh.boundary_config['type'] == 'dirichlet':
+ return dirichlet_boundary(projection=projection,
+ )
else:
raise Exception('Not implemented!')
return boundary
@@ -71,9 +98,9 @@ def enforce_boundary(num_ghost_cells=None):
def enforce_boxplot_boundaries(func):
def boxplot_boundary(self, *args, **kwargs):
bounds = np.array(func(self, *args, **kwargs))
- if self._mesh.boundary == 'periodic':
+ if self._mesh.boundary['type'] == 'periodic':
return tuple(periodic_boundary(projection=bounds,
- num_ghost_cells=1))
+ config=self._mesh.boundary_config))
else:
raise Exception('Not implemented!')
return boxplot_boundary
@@ -82,9 +109,8 @@ def enforce_boxplot_boundaries(func):
def enforce_fold_boundaries(func):
def fold_boundaries(self, *args, **kwargs):
folds = np.array(func(self, *args, **kwargs))
- if self._mesh.boundary == 'periodic':
+ if self._mesh.boundary['type'] == 'periodic':
return folds % self._mesh.num_cells
else:
raise Exception('Not implemented!')
return fold_boundaries
-
diff --git a/scripts/tcd/Mesh.py b/scripts/tcd/Mesh.py
index abd4baa90b96f6eb34bc7e4f0d37f8836911a826..2129342547296368b9a404da523ab27084d8e5dc 100644
--- a/scripts/tcd/Mesh.py
+++ b/scripts/tcd/Mesh.py
@@ -11,6 +11,9 @@ from typing import Tuple
from functools import cache
+BOUNDARY_TYPES = ['periodic', 'dirichlet']
+
+
class Mesh:
"""Class for mesh.
@@ -45,7 +48,7 @@ class Mesh:
"""
def __init__(self, num_cells: int, left_bound: float, right_bound: float,
- num_ghost_cells: int = 1) -> None:
+ boundary_config: dict, num_ghost_cells: int = 1) -> None:
"""Initialize Mesh.
Parameters
@@ -80,7 +83,18 @@ class Mesh:
'an exponential of 2')
self._left_bound = left_bound
self._right_bound = right_bound
- self.boundary = 'periodic'
+
+ self._boundary_config = boundary_config
+ boundary_type = self._boundary_config.get('type', 'periodic')
+ if boundary_type not in BOUNDARY_TYPES:
+ raise ValueError(f'The boundary type "{boundary_type}" '
+ f'is not implemented')
+ self._boundary_config['type'] = boundary_type
+
+ @property
+ def boundary_config(self) -> dict:
+ """Return boundary config."""
+ return self._boundary_config
@property
def num_cells(self) -> int:
@@ -127,7 +141,9 @@ class Mesh:
return {'num_cells': self._num_cells,
'left_bound': self._left_bound,
'right_bound': self._right_bound,
- 'num_ghost_cells': self._num_ghost_cells}
+ 'num_ghost_cells': self._num_ghost_cells,
+ 'boundary_config': self._boundary_config
+ }
def random_stencil(self, stencil_len: int) -> Mesh:
"""Return random stencil.