Skip to content
Snippets Groups Projects
Commit ab660bde authored by Laura Christine Kühle's avatar Laura Christine Kühle
Browse files

Added Dirichlet boundary.

parent 520bb6f2
No related branches found
No related tags found
No related merge requests found
......@@ -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
......
......@@ -22,7 +22,7 @@ class Basis(ABC):
Attributes
----------
polynomial degree : int
polynomial_degree : int
Polynomial degree.
basis : ndarray
Array of basis.
......
......@@ -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
......@@ -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.
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment